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The  electron-nuclear  dynamics  in  chemical  processes  is  described  by  a method 
founded  on  the  Time-Dependent  Variational  Principle  (TDVP).  In  addition,  to  avoid 
redundancies  in  the  parametrization  of  the  electronic  wave  function,  coherent  states  (CS) 
are  used.  The  equations  resulting  from  this  treatment  enable  us  to  describe  the  dynamics 
in  general  molecular  system,  since  it  does  not  presume  a priori  knowledge  of  the  potential 

energy  surface  of  the  molecular  system.  The  theory  is  called  Electron  Nuclear  Dynamics 
(END). 

In  order  to  make  this  time-dependent  method  based  upon  TDVP  and  CS  practical 
and  still  realistic,  we  treat  the  nuclei  in  the  limit  of  narrow  Gaussian  wavepackets  (which 
corresponds  to  classical  nuclei)  and  use  a single  determinant  to  describe  the  electrons. 
We  name  this  approach  END-SD-FGWP. 

We  use  the  END-SD-FGWP  approach  to  study  the  dynamics  of  charge  transfer  in 
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the  systems  Li^-H— Li  ^ Li— H-Li^  and  Li+-CN-Li  ^ Li-CN-Li-".  The  collision  of  a 
proton  with  H,  He,  and  H2  is  also  investigated  by  the  END-SD-FGWP  approach.  Several 
properties  of  these  collisions  are  compared  directly  with  the  experimental  results. 

Since  we  are  interested  in  the  application  of  the  END  formalism  to  large  systems  we 
propose  one  more  approximation,  namely,  the  neglect  of  the  differential  diatomic  over- 
lap (NDDO),  and  introduce  the  END-SD(NDDO)-FGWP  approach.  The  MNDO/AMl 
parametrization  is  used  for  the  NDDO  approach.  The  computational  savings  are  twofold: 
a)  the  number  of  integrals  computed  are  reduced  by  several  orders  of  magnitude  com- 
pared to  the  ab  initio  approach,  b)  the  equations  of  motion  are  simplified,  and  c)  the  high 
frequency  motions  of  the  core  electrons  is  removed  by  an  effective  core. 
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CHAPTER  1 
INTRODUCTION 


Time-dependent  methods  for  solving  the  Schrodinger  equation  were  introduced  in 
quantum  mechanics  since  its  beginning.  For  instance,  already  in  the  1930s  Dirac  sug- 
gested the  time-dependent  Hartree-Fock  method  [1].  However,  due  to  the  theoretical 
and  computational  complexity  of  time-dependent  methods,  their  development  and  ap- 
plication to  molecular  systems  have  been  delayed  by  several  decades  in  comparison  to 
time-independent  methods.  Despite  these  difficulties,  time-dependent  methods  for  molec- 
ular quantum  dynamics  have  gained  in  importance  in  the  last  two  decades.  Certainly, 
the  interest  in  time-dependent  solutions  of  the  Schrddinger  equation  is  very  much  due  to 
the  development  in  recent  years  of  novel  experimental  techniques  that  make  it  possible 
to  measure  quantum  state-specific  properties  of  rather  complicated  systems  and  to  probe 
the  dynamics  of  chemical  reactions  at  the  molecular  level  for  very  short  times.  However, 
a solid  theoretical  understanding  of  non-stationary  processes  is  highly  desirable,  and  due 
to  the  development  of  new  computational  techniques  and  powerful  computer  hardware, 
assumptions,  hypotheses  and  theories  describing  these  processes  can  now  be  explored. 

The  time-dependent  approach  to  quantum  dynamics  offers  a number  of  advantages 
over  time-independent  methods.  For  example,  in  a reactive  scattering  process,  plots  of 
the  particle  density  as  it  evolves  throughout  the  collision  are  easily  interpreted  using 
simple  concepts  of  the  underlying  physics.  These  plots  can  be  used  to  improve  one’s 
understanding  of  mechanisms  that  are  important  during  a reaction  or  in  other  physical- 
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chemical  processes.  From  a theoretical  point  of  view,  the  time-dependent  approach  makes 
it  straightforward  to  include  the  effect  of  intense  laser  fields  or  external  time-dependent 
heat  baths,  where  the  dynamics  is  treated  approximately  by  classical  mechanics,  statistical 
or  stochastic  methods.  It  is  also  well  adapted  to  improvement  in  the  dynamics  for  those 
degrees  of  freedom  that  are  described  classically,  by  the  wavepacket  propagation,  as 
well  as  those  given  by  explicit  quantum  mechanical  treatments.  As  a result,  the  END 
and  other  time-dependent  approaches  offer  great  flexibility  for  the  choice  of  various 
approximations  and  is  also  capable  of  a dynamically  accurate  description.  It  is  also 
possible  to  use  within  these  approaches  some  of  the  experience  acquired  for  decades 
with  time-independent  methods  to  provide  a practical  and  realistic  description  of  the 
chemical  system  while  at  the  same  time  keeping  an  accurate  description  of  the  dynamics. 
It  also  allows  the  blending  of  classical,  wavepacket,  quantum  operator,  state  expansion 
and  numerical  techniques  necessary  for  solving  complex  dynamical  problems. 

We  present  a time-dependent  method  for  solving  the  Schrddinger  equation  and  explore 
its  capabilities  in  describing  charge  exchange  collisions,  molecular  vibrations,  and  electron 
transfer  dynamics.  The  method  presented  here  attempts  to  describe  the  simultaneous 
electron  and  nuclear  dynamics.  It  has  been  developed  during  the  last  half  decade  and  has 
the  very  nice  feature  of  not  requiring  potential  energy  surfaces  of  the  molecular  system. 
It  allows  both  electrons  and  nuclei  to  be  treated  quantum  mechanically.  The  dynamical 
equations  of  motion  are  obtained  using  the  time-dependent  variational  principle  (TDVP) 
and  the  wavefunctions  are  parametrized  using  coherent  states  (CS).  The  combination  of 
TDVP/CS  gives  rise  to  a new  time-dependent  approach  to  molecular  dynamics  called 
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Electron  Nuclear  Dynamics  (END).  The  END  formalism  provides  dynamical  equations 
that  look  simple  and  are  easy  to  interpret.  The  computer  implementation  of  this  approach 
is  not  trivial,  and  in  order  to  make  it  practical  some  approximations  are  necessary.  Based 
upon  the  experience  acquired  with  time-independent  methods,  the  nuclei  are  treated 
classically  or,  more  appropriately,  in  the  limit  of  narrow  Gaussian  wavepackets,  and 
the  electrons  are  described  by  a single  determinant.  The  theory  for  quantum  nuclei  and 
multi-determinantal  description  of  the  electrons  has  already  been  developed  [2],  but  its 
implementation  is  not  yet  done.  So  we  think  it  would  be  interesting  and  fruitful  to  explore 
in  some  detail  the  single  determinant/classical  nuclei  model.  Certainly,  this  simplest 
model  for  the  END  formalism  may  be  inadequate  for  describing  some  problems,  mainly 
the  ones  involving  many  product  channels  (multi-reference  character)  and  tunneling 
dominant  processes  [3].  However,  since  the  END  formalism  does  not  require  the 
knowledge  of  the  potential  energy  surface(s)  of  the  molecular  system  it  can  be  applied  to 
a large  variety  of  systems  and  dynamical  problems  in  an  unified  manner.  For  instance, 
it  can  be  used  to  study  scattering  problems  (total,  differential  and  state  to  state  reactive 
cross  sections),  charge  transfer  dynamics  (rates  and  mechanisms),  molecular  vibrations 
(interpretation  of  transition  state  spectroscopy  experiments  and  prediction  of  reaction 
rates).  It  permits  us  to  follow  chemical  reactions,  predict  and  interpret  Raman  spectra,  and 
many  other  dynamical  processes.  Future  implementation  of  time-dependent  field  in  the 
molecular  Hamiltonian  would  make  the  END  formalism  able  to  describe  photodissociation 
processes,  laser  chemistry,  dynamical  nuclear  magnetic  resonance  (NMR),  etc.  It  should 
be  noted  however,  that  a fully  ab  initio  treatment  of  large  molecular  systems  using 
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the  END/single  determinant/classical  nuclei  approach  is  not  yet  practical.  Therefore,  in 
addition  to  exploring  the  ab  initio  version  for  describing  reactive  collisions,  molecular 
vibrations,  and  electron  transfer  in  small  systems,  we  propose  to  employ  the  neglect  of 
the  differential  diatomic  overlap  (NDDO)  approximation  [4].  It  is  our  expectation  that 
the  NDDO  approximation  on  top  of  the  END/single  determinant/classical  nuclei  model 
will  decrease  the  computational  demand  in  many  ways. 

a)  The  number  of  molecular  integrals  and  their  gradients  involved  would  be  several  orders 
of  magnitude  smaller  than  the  ab  initio  version. 

b)  Since  it  is  commonplace  to  use  effective  potential  for  the  core  electrons  it  should,  in 
addition  to  decreasing  the  number  of  degree  of  freedom  (differential  equations  needed 
to  be  solved)  also  eliminate  the  high  frequency  motions  of  the  core  electrons.  This  will 
allow  the  integrator  of  the  differential  equation  to  take  larger  steps  decreasing  by  several 
orders  of  magnitude  the  time  to  solve  the  differential  equations. 

c)  In  general,  minimum  basis  sets  are  the  working  frame  for  the  NDDO  approximation, 
which  decreases  significantly  the  number  of  equations  of  motion. 

The  outline  of  this  work  is  as  follows.  Chapter  2 describes  some  of  the  available 
time-dependent  methods  for  molecular  dynamics.  The  END  formalism  is  introduced  in 
Chapter  3 with  the  NDDO  approximation  and  concluded  with  the  END-NDDO  equations. 
In  Chapter  4 we  present  the  results,  analysis  and  discussion  of  the  reactive  scattering  of 
on  H,  He,  and  H2.  Chapter  5 is  devoted  to  the  study  of  the  intramolecular  charge  transfer 

systems  y-CN-Li  ^ Li-CN-Li-"  and  y-H— Li  ^ Li— H-Li-".  The  conclusions  are 
presented  in  Chapter  6. 


CHAPTER  2 

TIME-DEPENDENT  METHODS 


The  existing  time-dependent  methods  that  are  relevant  to  our  work  are  described  in 
this  chapter.  The  aim  is  to  provide  a general  view  of  these  methods  with  emphasis  on 
performance,  limitations,  and  comparisons  among  them.  It  seems  that  the  best  way  to  give 
an  overview  of  these  time-dependent  methods  is  to  provide  some  kind  of  classification 
that  compare  their  main  assumptions.  In  order  to  provide  such  a classification,  the  concept 
of  potential  energy  (hyper)surface  (PES)  is  fundamental.  For  practical  purposes  a PES 
is  a function  that  gives  the  value  of  the  electronic  energy  of  a group  of  N interacting 
atoms  with  respect  to  their  37V  — 6 relative  position  coordinates.  For  each  electronic 
state  of  the  molecular  system  we  have  one  PES  which  contains  information  about 
the  isolated  reactant  and  product  species,  their  long-range  interactions,  their  molecular 
deformations  which  lead  to  the  breaking  and  forming  of  chemical  bonds,  the  geometries 
and  properties  of  transient  reaction  intermediates,  and  energy  barriers  for  going  from  one 
of  those  intermediates  to  another.  As  we  can  see,  the  task  involved  in  obtaining  such 
a (hyper)surface  for  polyatomic  molecules  is  tremendous.  Highly  correlated  electronic 
structure  calculations  with  very  large  basis  sets  are  necessary  in  order  to  provide  an 
accurate  description  of  the  different  regions  of  the  PES.  Once  the  PES  has  been  mapped, 
a fitting  procedure  will  be  employed  to  give  the  functional  form  of  the  PES.  This  fitting 
procedure  is  very  time  consuming  and  tedious,  and  it  is  still  more  an  art  than  a scientific 
procedure.  The  functional  form  of  a PES  can  be  quite  complicated.  As  evidence  for  that 
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we  can  mention  that  even  for  a triatomic  molecule  like  H3,  an  attempt  to  use  the  well 
developed  algebra  program  (REDUCE)  to  evaluate  the  first,  second,  third,  and  fourth 
derivatives  analytically  failed  due  to  the  algebraic  complexity  [5].  We  classify  the  time- 
dependent  methods  into  two  categories  a)  PES  based  methods,  and  b)  PES  free  methods. 

PES  Based  Time-Dependent  Methods 

Once  the  potential  energy  surface  for  the  molecular  system  has  been  determined, 
there  are  basically  three  approaches  to  perform  dynamics  on  this  surface,  namely,  a) 
the  classical  trajectory  method,  b)  the  semiclassical  method,  and  c)  the  quantum  time- 
dependent  method.  The  given  order  a),  b),  and  c)  is  related  to  the  increasing  theoretical 
and  implementation  difficulties,  computational  demand,  and  accuracy. 

The  widely  used  classical  trajectory  approach  [6-8]  employs  Newton’s  equation 
of  motion  to  treat  the  dynamics  of  the  nuclei  on  the  PES.  This  approach  is  easily 
implemented  and  the  interpretation  of  the  results  are  simple.  However,  the  computational 
demand  can  be  large  for  polyatomic  systems,  and  a good  statistical  analysis  program  is 
necessary.  A detailed  classical  trajectory  simulation  requires  an  extensive  sampling  in 
the  phase  space  (50  000  to  100  000  trajectories  for  a triatomic  system)  for  each  initial 
condition  (collision  energy,  vibrational  and  rotational  initial  states).  This  computational 
demand  is  not  the  only  drawback  of  the  classical  method.  It  neglects  quantum  effects  like 
tunneling,  interference,  penetration  (generalized  tunneling),  zero-point  motion,  resonance 
behavior,  etc,  which  can  be  fundamental  for  chemical  encounters.  These  deficiencies 
represent  an  inherent  obstacle  for  the  classical  trajectory  approach  to  provide  qualitative 
or  quantitative  results  for  chemical  processes.  Despite  these  problems,  the  classical 
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trajectory  approach  has  been  applied  to  a widely  variety  of  dynamical  problems  including 
scattering  of  atoms  or  molecules  by  solid  surfaces  [9]  and  bimolecular  reactions  [8],  In 
most  of  these  applications  (semi)empirical  potential  energy  surfaces  are  used,  such  as 
DIM  (Diatomic-in-Molecules),  MIM  (Molecule-in-Molecules),  LEPS  (London-Eyring- 
Polany-Sato),  etc  [10,  11]. 

Semiclassical  methods  have  been  designed  to  overcome  some  of  the  flaws  in  the 
classical  description.  They  assume  that  the  quantum  effects  can  be  treated  as  corrections 
to  the  classical  motion,  based  on  the  fact  that  the  nuclei  are  relatively  heavy  compared 
to  the  electronic  mass.  The  semiclassical  approaches  have,  in  principle,  the  advantage 
of  incorporating  some  of  the  important  quantum  effects  while  maintaining  the  simple 
implementation.  However,  the  generalization  of  the  semiclassical  methods  is  not  simple, 
and,  as  result,  each  system  requires  a special  implementation,  and  the  validity  of  the 
approximations  is  often  not  clear.  The  perhaps  most  developed  semiclassical  dynamical 
method  has  been  proposed  by  Heller  [12—15]  and  employs  the  Gaussian  wavepacket 
approximation.  That  is,  the  quantum  wavepacket  is  constrained  to  be  of  a Gaussian 
form.  This  semiclassical  method  has  been  widely  used  in  many  applications  with  some 
success  [12-15].  Another  well-developed  semiclassical  method  that  uses  the  so  called 
eikonal  approach  to  dynamics  has  been  proposed  by  Micha  [16].  Comparisons  of  the 
results  provided  by  these  two  semiclassical  methods,  namely  the  Gaussian  wavepacket 
and  the  eikonal  approach,  have  shown  that  they  differ  only  in  some  minor  details  [17,  18]. 
These  semiclassical  methods  have  been  applied  to  describe  photodissociation  processes, 
reactive  scattering  and  bimolecular  reactions  [15]. 
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Among  the  various  quantum  time-dependent  methods  on  PES  two  have  been  widely 
used,  namely,  the  ’exact’  time-dependent  quantum  mechanical  and  the  quantum  time- 
dependent  self-consistent  field  (TDSCF)  methods  [19,  20].  Even  the  so  called  ’exact’ 
time-dependent  methods  do  use  some  approximations,  which  lead  to  errors  that  can,  in 
principle,  be  made  as  small  as  desired.  The  ’exact’  time-dependent  methods  can  be 
classified  according  to  the  representation  of  the  wavefunction  and  the  algorithm  for 
the  time  evolution.  We  distinguish  two  approaches  to  wavefunction  representation, 
namely,  discretization  on  a grid  of  points,  sometimes  referred  to  as  DVR  (discrete 
variable  representation),  and  expansion  in  a basis  set.  In  general,  the  numerical  cost 
of  the  discretization  and  basis  set  expansion  representations  scales  as  0{V^)  and  0{N^), 
respectively,  where  V is  the  volume  of  phase  space  containing  the  system  and  N is  the 
number  of  expansion  functions.  There  are  basically  four  algorithms  in  use  for  the  time 
evolution:  SOD  (second-order  differencing),  SO  (split  operator),  SIL  (short-time  iterative 
Lanczos),  and  CE  (Chebychev  expansion).  For  a detailed  presentation  and  comparisons  of 
the  different  ’exact’  time-dependent  methods,  see  references  [20,  21].  Applications  of  the 
’exact’  time-dependent  methods  have  been  restricted  to  at  most  three  degrees  of  freedom 
and  include  mainly  photodissociation  processes,  absorption  and  photodetachment  spectra, 
reactive  scattering  (atom-diatom),  and  unimolecular  reactions  [22].  The  TDSCF  method 
assumes  variables  separation  such  that  each  degree  of  freedom  is  treated  in  the  mean  field 
of  the  remaining  degrees  of  freedom.  In  the  static  SCF  this  approximation,  also  called 
the  independent  particle  model,  leads  to  the  lack  of  correlation  between  the  degrees  of 
freedom  of  the  system.  However,  in  the  case  of  TDSCF  method  energy  is  allowed  to  flow 
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between  degrees  of  freedom  leading  to  some  correlation  among  them.  As  a result,  the 
TDSCF  scheme  is  expected  to  be  accurate  for  calculations  of  such  averaged  properties  as 
lifetimes,  mean  fragment  energies,  energy  flows,  etc.  In  fact,  the  dissociation  dynamics 
of  the  l2He  van  der  Waals  complex  has  been  examined  by  ’exact’  and  TDSCF  methods. 
It  has  been  shown  that  for  this  system  the  TDSCF  scheme  gives  results  in  excellent 
agreement  with  ’exact’  method  when  the  calculated  properties  are  obtained  during  a 
period  of  time  much  larger  than  the  vibrational  periods  and  so  long  as  the  excitation 
level  is  not  too  high.  The  TDSCF  approach  also  reproduces  well  the  dephasing  events 
in  the  dynamics  and  gives  very  good  results  for  the  autocorrelation  function  [19,  20]. 

We  can  summarize  this  section  by  saying  that  the  greatest  drawback  of  these  methods 
(classical,  semiclassical,  ’exact’,  TDSCF)  is  the  need  for  one  or  more  potential  energy 
surfaces.  In  the  case  of  classical  methods  the  lack  of  any  quantum  effects  can  cause 
problems,  the  semiclassical  approaches  are  too  system  specific,  and  the  ’exact’  schemes 
are  very  numerically  intensive  and  are  restricted  to  problems  with  few  degrees  of  freedom. 
No  doubt  these  methods  have  their  advantages  and  their  contribution  to  the  understanding 
of  time-dependent  processes  must  be  emphasized. 

PES  Free  Time-Dependent  Methods 


We  focus  in  this  section  on  time-dependent  methods  that  do  not  require  potential 
energy  surfaces.  Among  these  methods  we  intend  to  discuss  the  time-dependent  methods 
that  fall  into  these  three  classes:  a)  close-coupling,  b)  molecular  dynamics,  and  c)  time- 
dependent  Hartree-Fock. 
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Close-coupling  Methods 

The  close-coupling  methods  [23-25]  consist  in  projecting  the  time-dependent 
Schrodinger  equation  on  a basis.  The  electronic  states  are  described  by  a Cl  (config- 
uration interaction)  expansion  where  the  time-dependence  is  put  into  the  Cl  coefficients. 
The  nuclear  dynamics  are  not  treated  explicitly  and  are  constrained  to  follow  a priori 
trajectories,  such  as  linear  or  Coulombic  trajectories. 

The  description  of  the  electronic  states  by  a Cl  wavefunction  imposes  a heavy 
computational  demand  restricting  the  application  of  the  close-coupling  methods  to  simple 
one  electron  (or  pseudo  one-electron)  systems.  Recent  extensions  to  two-electron  systems 
seems  not  to  be  very  practical  [26]. 

The  prescribed  trajectories  for  the  nuclei  limit  the  application  of  the  close-coupling 
methods  to  high  energy  (above  1.0  keV)  collisions.  It  also  makes  the  application  of 
these  methods  for  the  calculation  of  differential  cross  sections  other  than  the  smallest 
scattering  angles  quite  doubtful. 

Molecular  Dynamics  Methods 

Molecular  dynamics  (MD)  methods  assume  that  the  atomic  dynamics  is  governed 
by  the  laws  of  classical  mechanics.  That  is,  MD  methods  neglect  quantum  effects  on 
the  nuclear  motion  and  the  electrons  are  considered  as  being  always  in  the  ground- 
state  corresponding  to  the  instantaneous  nuclear  configuration.  As  a result  the  explicit 
knowledge  of  the  interatomic  potential  is  required.  Obtaining  this  interatomic  potential  is 
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a difficult  many-body  problem  and  different  ways  of  calculating  this  potential  distinguish 
the  different  MD  methods. 

The  most  widely  used  MD  methods  are  the  ones  that  employ  (semi)empirical  force 
fields  for  describing  the  interatomic  potential  [27].  These  MD/force-field  methods  can 
be  combined  with  perturbation  methods  to  provide  thermodynamical  properties  [28],  in 
addition  to  structural  (geometric  and  conformational)  information.  They  can  be  applied 
to  a great  variety  of  systems  including  large  biological  molecules  [29]. 

Some  of  the  other  MD  formulations  evolve  computing  the  interatomic  potential 
by  solving  the  Schrddinger  equation  in  some  approximate  manner  (e.g.  ab  initio  or 
semiempirical  Hartree-Fock  (HF),  density  functional  (DF)  theory,  etc)  for  each  nuclear 
configuration  during  the  time  evolution  [30,  31]. 

Time-Dependent  Hartree-Fock  Methods 

Despite  the  fact  that  the  time-dependent  Hartree-Fock  (TDHF)  method  was  suggested 
by  Dirac  already  in  1930s  [1]  only  recently  has  this  method  been  developed  and  applied 
to  reactive  and  inelastic  scattering  [32-35]  and  molecular  problems  [36,  37]. 

In  TDHF  the  quantum  dynamics  of  the  electrons  represented  by  a single  Slater 
determinant  is  coupled  to  a classical  description  of  the  nuclear  dynamics  in  a self- 
consistent  fashion.  The  differences  between  some  of  the  TDHF  implementations  are 
related  to  the  way  the  Hartree-Fock  equation  is  approximately  solved  and  to  the  way  that 
the  translation  of  the  basis  set  with  the  nuclei  is  handled.  This  is  related  to  the  electron 
translation  factor  (ETF)  problem. 
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The  NDDO  approximation  has  been  used  to  simplify  the  TDHF  equations  and  applied 
in  the  calculation  of  molecular  properties,  such  as  atomic  charges,  structure,  dipole 
moments,  etc,  during  time  evolution  of  LiH,  H2O,  and  CH2O  molecules  [36].  The  density 
functional  theory  (DFT)  using  the  local  density  approximation  (LDA)  has  been  employed 
to  compute  the  Hartree-Fock  potential  for  the  TDHF  method.  This  approach  has  been 
used  in  the  dynamics  of  Na4  cluster  for  the  singlet  and  triplet  electronic  states  [37].  Both 
approaches,  namely,  TDHF-NDDO  and  TDHF-DFT,  do  not  consider  the  effects  of  the 
translation  of  the  basis  set  with  the  nuclear  motion  upon  the  electron  dynamics.  This 
means  that  the  electron  translation  factors  are  not  properly  taken  into  account  by  these 
approaches,  neither  are  the  electron-nuclear  couplings. 

The  proper  derivation  of  TDHF  equations  with  the  inclusion  of  the  electron  translation 
factors  has  been  done  using  the  density  matrix  instead  of  the  orbital  coefficients  as  the 
time -dependent  quantities  [34].  However,  no  applications  of  this  method  have  yet  been 
made,  and  only  after  the  neglect  of  the  electron  translation  factors  (RTF’s)  was  this 
approach  applied  to  the  H'*’  + H collision  [34].  It  is  expected  that  in  a short  time 
applications  including  the  ETF’s  will  appear  and  we  shall  be  able  to  see  the  magnitude 
of  the  errors  caused  by  the  neglect  of  these  effects. 

In  summary: 

1 . PES  based  time-dependent  methods  have  the  drawback  of  requiring  a priori  knowl- 
edge of  the  interparticle  potentials.  This  requirement  limits  the  applicability  of  these 
methods  to  small  systems  (just  a few  degrees  of  freedom)  and  low  energy  processes 
(high  energy  implies  in  a large  number  of  excited  PES’s), 
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2.  PES  free  time-dependent  methods  are  very  promising,  but  general  computer  codes, 
extensive  testing  of  the  approximations,  and  a solid  foundation  of  the  formal  devel- 
opments are  still  required. 

As  a conclusion  to  this  chapter  we  can  say  that  there  are  still  a lot  of  formal 
and  numerical  investigations  to  be  performed  on  time-dependent  methods  for  molecular 
quantum  dynamics.  In  the  next  few  chapters  we  shall  be  able  to  contribute  to  the 
investigation  of  a new  time-dependent  method. 


CHAPTER  3 

THE  THEORETICAL  APPROACH 


In  the  last  chapter  we  expressed  our  concern  for  potential  energy  surface  (PES)  based 
time-dependent  methods  being  limited  to  systems  where  these  PES’s  are  known.  In  this 
chapter  we  present  the  equations  for  a recently  proposed  time-dependent  approach  to 
molecular  dynamics  called  Electron  Nuclear  Dynamics  (END). 

The  time-dependent  variational  principal  (TDVP)  is  used  to  obtain  an  approximate 
time  evolution  of  the  molecular  system.  In  order  to  obtain  the  equations  of  motion  it  is 
necessary  to  decide  how  the  electrons  and  nuclei  will  be  treated.  From  the  experience 
with  time-independent  methods  and  some  of  the  time-dependent  methods,  like  TDSCF 
and  TDHF,  it  seems  that  treating  the  electrons  at  the  single  determinant  level  with  classical 
nuclei  should  provide  a very  realistic  approach  for  a wide  variety  of  chemical  processes. 
Since  the  classical  nuclei  approximation  can  be  obtained  as  the  narrow  width  limit  of 
a frozen  Gaussian  wave  packet  treatment,  we  denote  the  method  resulting  from  these 
approximations  as  END-SD-FGWP.  That  is,  END  stands  for  Electron  Nuclear  Dynamics 
from  the  approximate  dynamical  equations  obtained  from  the  TDVP,  SD  is  the  Single 
Determinantal  description  of  the  electrons,  and  FGWP  is  the  classical,  or  narrow  width 
limit  of  Frozen  Gaussian  Wave  Packet  representation  of  the  nuclei.  The  coherent  state 
(CS)  for  the  electrons  is  the  Thouless  single  determinant  representation,  which  gives  a 
suitable  parameterization. 

In  order  to  make  the  END-SD-FGWP  approach  suitable  for  studying  dynamics  of 
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large  molecular  systems,  we  employ  the  neglect  of  diatomic  differential  overlap  (NDDO) 
approximation. 

We  have  divided  this  chapter  into  five  sections,  namely,  the  END  formalism,  the 
END-SD-FGWP  method,  the  NDDO  approximation,  the  END-SD(NDDO)-FGWP  model, 
and  implementation.  In  the  first  section  the  CS,  the  TDVP,  and  the  dynamical  equations 
are  introduced.  In  the  NDDO  section  we  have  a subsection  where  we  present  the  Austin 
model  1 (AMI)  realization  of  the  NDDO  approximation. 

The  END  Formalism 

We  present  a formalism  that  treats  the  electron  nuclear  dynamics  without  the  need  of 
a priori  knowledge  of  potential  energy  surfaces  or  interatomic  potentials.  In  other  words, 
in  this  formalism  only  the  Hamiltonian  (e.g.,  Coulombic)  and  the  initial  conditions  (initial 
nuclear  positions  and  momenta,  for  instance)  are  needed  for  the  calculation  to  proceed. 
This  formalism  relies  upon  the  time-dependent  variational  principal  (TDVP)  of  a quantum 
state  describing  the  system  through  a group  theoretical  coherent  state  (CS). 

In  the  next  two  subsections  we  present  the  dynamical  equations  for  a molecular 
system  derived  using  the  TDVP  and  the  CS  representation.  The  actual  derivation  has 
been  presented  elsewhere  [38,  39]  and  will  not  be  repeated  here. 

Coherent  States 

We  now  present  the  parametrization  of  the  wavefunction  in  terms  of  coherent  states 
(CS).  A coherent  state  [38,  40]  of  a molecular  wave  function  |C)  is  a continuous  function 
of  the  labeled  set  of  electronic  and  nuclear  parameters  (.  This  set  of  parameters  has  a 


16 


positive  measure  d(  such  that  the  resolution  of  the  identity 

j ICHCWC  = 1 


(3.1) 


holds. 

This  choice  of  representing  the  wave  function  is  very  important  since  it  removes  any 
parameter  redundancies  and  leads  to  a simple  interpretation  of  the  dynamical  equations. 
The  removal  of  the  redundant  parameters  is  essential  to  avoid  complications  with  the 
phase  space  that  could  lead  to  artificial  singularities  and  difficulties  with  integration.  A 
coherent  state  representation  also  leads  to  a natural  division  of  the  parameters  into  (gen- 
eralized) canonical  coordinates  and  conjugate  momenta.  Thus,  they  form  a generalized 
phase  space  with  the  expectation  value  of  the  energy  being  the  Hamiltonian  function  for 
the  time  evolution  of  the  parameters.  As  a result,  the  equations  for  the  CS  parameters 
form  a classical  Hamiltonian  system  facilitating  their  physical  interpretation. 

The  generalization  to  a multi-determinantal  description  would  follow  the  same  general 
approach  and  has  been  done  elsewhere  [2] . We  are  interested  in  those  coherent  states  that 
are  related  to  compact  Lie  groups  and  to  a single  determinant.  Using  the  group  theoretical 
approach  we  can  eliminate  the  redundancy  of  the  single  determinant  representation.  Let 
us  assume  that  we  have  an  orthonormal  spin  orbital  basis  set  made  up  of  K vectors 

and  N particles  which  will  be  described  by  a single  determinant.  This 
determinant  can  be  found  as  a unitary  transformation  acting  on  the  orbital  expansion 
coefficients  of  some  reference  single  determinant  |^o)-  We  define  this  reference  state 
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l^o)  = n^^K^)  (3.2) 

/=1 

where  b]  (6,)  are  the  fermion  creation  (annihilation)  operators  of  the  spin  orbital  basis 
used.  The  reference  state  is  a lowest  weight  state  for  the  irreducible  representation 
[l'^0(^“^)]  of  the  group  U{K)  with  generators  bjbj  because  the  weight  or  number 
operators  b]b{  have  eigenvalues  1 for  * = 1, . . . , iV  and  0 for  i = iV  + 1, . . . , ii'.  As  we 
can  see  we  have  divided  our  set  of  K spin  orbitals  into  two  subsets  N and  K — N.  We 
now  denote  any  quantity  q as: 

1.  q if  q refers  to  the  set  of  all  K spin  orbitals; 

2.  q*  if  q refers  to  the  subset  of  N spin  orbitals;  and 

3.  9°  if  q refers  to  the  subset  of  K - N spin  orbitals. 


In  the  context  of  a determinantal  state  this  division  of  the  set  of  K spin  orbitals  into 
the  two  subsets:  N and  K — N takes  the  meaning  of  occupation,  where  the  subset  N 
is  called  occupied  and  K — N is  the  unoccupied  or  virtual  subset.  In  the  general  case, 
however,  these  subsets  N and  K — N will  be  bullet  and  open  circle,  respectively.  Also, 
the  following  notation  will  be  used  when  dealing  with  matrices: 

/ “11  •••  ^IN  “l(7V+l)  •••  ^'iK  \ 


u = 


u 


^N1 

(AT+1)1 


u 


n 

K\ 


u 


u 

n 


NN 


u 


It 

KN 


u 


N{N^\) 


(TV-fl)A^  ^(Ar-hi)(AT+i) 


u 


A'(iV+l) 


u 


NK 


u 


u 


KK 


(3.3) 


^NxN 

If" 

^K-NxN 


^NxK-N  \ _ 
^K-NxK-N 


• • • 


/ 
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We  now  assume  that  there  is  an  t/  matrix  which  transforms  the  spin  orbital  basis 
to  the  occupied  and  virtual  spin  orbitals,  such  that,  all  redundancies  due  to  unitary 
transformations  among  the  occupied  (virtual)  spin  orbitals  are  removed,  that  is. 


(3.4) 


Then,  the  coherent  state  representation  of  a determinantal  state  takes  the  form: 


N 


1=1 


= n ( E + E 1 1'"“^) 

i=l  \i=l  ji=7V+l 


N / N 


K N 


= n I E 1 + E E f'yikVtr'  I f'l*  I \o“c) 

j=N+l  k=l 


i=l  \/=l 


(3.5) 


“HK+  E 

i=l  V >=A^+lib=l 


=a 


n('+  E 

i=i  \ y=jv+iJb=i  / 1=1 
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and  since  we  work  with  an  unnormalized  coherent  state  we  have  that 


N K 


=n  n 


t=l  j=N-\-l 


(3.6) 


N K 


1=1  j=N+l 


where  ^ 


The  unnormalized  determinantal  wavefunction  can  be  expressed  in  terms  of  the 
orthonormal  spin  orbital  basis  set  and  the  dynamical  variables  {zji,Zji} 

in  the  following  way 


where  the  dynamical  spin  orbitals 

K 

Xi  = V’i  + (1  < * < N)  (3.8) 

j=N+i 

are  nonorthogonal,  but  linearly  independent.  The  corresponding  unoccupied  dynamical 
orbitals  may  be  chosen  as 


and  although  mutually  nonorthogonal  they  are  orthogonal  to  the  occupied  space. 


\z)  = det{x,} 


(3.7) 


N 


(3.9) 
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The  Time-Dependent  Variational  Principle 


We  now  give  an  outline  of  how  the  time-dependent  variational  principle  (TDVP)  can 
be  used  to  provide  us  with  an  approximate  time  evolution  of  the  quantum  mechanical 
system.  There  are  several  forms  of  the  TDVP  [1,  41,  42],  but  they  are  all  equivalent 
when  the  trial  wave  function  is  described  by  complex  parameters  and  is  analytic  in  this 
parameter  space  [43]. 

The  TDVP  approach  consists  in  minimizing  the  variations  of  the  quantum  mechanical 
action  functional  A,  that  is. 


f[(Cl|(IO) 


dt 


«ci)ic)  - (ci^io  Mcicr' = 0 


(3.10) 


where  \Q  is  quantum  state  of  the  system  labeled  with  a set  of  electronic  and  nuclear 
parameters  ( and  H is  the  Hamiltonian  of  the  system.  The  TDVP  approach  leads  then 
to  the  following  equation  for  the  time-dependent  state  |C) 


(3.11) 


This  is  the  time -dependent  Schrodinger  equation,  provided  that  \6()  is  allowed  to  vary 
over  all  Hilbert  space  and  the  right-hand  side  vanishes.  It  can  be  shown  [38]  the  right- 
hand  side  of  equation  (3.11)  vanishes  when  we  consider  explicitly  the  overall  phase 
factor  of  the  time-dependent  state  |(). 
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The  dynamical  equations  resulting  from  this  approach  are  [38,  39] 


E = 

0 


dE 

d(o 


(3.12) 


which  in  matrix  form  becomes 


C 0 
0 -C* 


c 

c* 


/ dE 

( w 

\ M 

V ac 


(3.13) 


and  the  elements  of  the  metric  matrix  defined  as 

^ a^InS 

^010  — p,. 


(3.14) 


dQ9C0 

where  the  energy  of  the  system  is  given  by  the  expectation  value  of  the  Hamiltonian 


H as  follows 


E{C\0  = 


{C\H\0 

(CIO 


(3.15) 


and  the  norm  is 


5(C*,C)  = (CIC). 


(3.16) 


It  also  should  be  mentioned  that  in  this  formulation  the  global  phase  factor  exp  (^7)  of 
the  wavefunction 


l^')  = exp  (n)  10 


(3.17) 


has  been  factorized  from  the  evolving  state  and  can  be  obtained  by  a simple  quadrature 
integration  of  the  following  differential  equation 


7 4E(<. 

a 


a,n5(c,o_^.ai^ 


\ CK 


dCa 


d\nSiC,0 

dCa 


da 


(3.18) 


-^(C*,C) 
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during  the  evolution.  This  global  phase  should  be  important  in  applications  where  the 
time-correlation  functions  need  to  be  computed. 

The  END-SD-FGWP  Method 

We  have  so  far  obtained  dynamical  equations  that  describe  the  dynamics  of  a general 
quantum  molecular  system.  As  mentioned  before,  in  order  to  obtain  equations  of  motion 
for  the  electron  nuclear  dynamics  that  are  practical  and  realistic,  some  approximations 
are  needed.  We  choose  a model  for  which 

(1)  the  electrons  are  described  by  a single  Slater  determinant,  and 

(2)  since  we  would  like  to  make  use  of  the  available  computer  codes  for  molecular 
integrals  we  use  a truncated  spin  orbital  basis  without  electron  translation  factors  (ETF’s). 
In  other  words,  the  basis  set  is  independent  of  the  nuclear  velocity  (or  momentum),  and 

(3)  use  a classical  description  of  the  nuclei.  In  mathematical  terms,  this  means  that  when 
employing  the  TDVP  approach  the  limit  of  narrow  Gaussian  wave  packets  is  taken.  As 
a result,  the  nuclear  parameters  become  the  positions  and  momenta  of  the  nuclei. 

TTiese  approximations  allow  us  to  obtain  the  equations  of  motion  for  END-SD-FGWP 
method. 

The  Equations  of  Motion 

Taking  the  limit  of  narrow  Gaussian  wave  packets  or  the  classical  description  of  the 
nuclei  permit  us  to  define  the  dynamic  variable  for  the  nuclei  as  [38,  39] 


Zk  = Rk  + iPk 


(3.19) 
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where  R and  P are  the  nuclear  coordinates  and  momenta,  respectively.  Since  we  have 
assumed  that  the  dependence  of  the  basis  on  the  nuclear  parameters  is  only  through  the 
nuclear  positions,  the  equations  of  motion  for  the  electron  and  nuclei  become  [38,  39] 


/ iC 

0 

iCR 

0 ^ 

(dEldz*\ 

0 

j 

-iC* 

0 

• ♦ 

dE/dz 

iC  R 

-iCl 

Crr 

-I 

k 

dEjdR 

0 

0 

I 

0 ) 

\p) 

\ dEjdP  J 

where 


d‘^\nS{z\R\P\z,R,P) 

dz*dXk 


R'=R,P'=P 


and 


CxkYi  = -2Im 


d^\nS{z\R',P\z,R,P) 

dX'dVi 


R'=R,P'z=P 


with  X and  Y standing  for  R or  P,  with  the  norm  S being 


(3.20) 


(3.21) 


(3.22) 


5 = S{z'\R',P',z,R,P)  = {z'\R!,P'\z,R,P) 
= det(/*  + z'^z) 


(3.23) 


In  order  to  implement  the  END-SD-FGWP  formalism  developed  so  far  we  need  to 
specify  how  to  construct  the  spin  orbital  basis.  Usually,  in  electronic  structure  calculation, 
a set  of  K localized  atomic  orbitals  are  used  to  define  the  spin  orbital  basis.  These  atomic 
orbitals  are  given  as  linear  contractions  of  Slater  (STO)  or  Gaussian  (GTO)  type  orbitals 
and  are,  in  general,  centered  on  the  nuclei.  These  orbitals  are  also,  in  general,  non- 
orthonormal. 


24 


Let  VL  be  a matrix  that  transforms  the  non-orthonormal  localized  atomic  orbital  basis 


to  the  orthonormal  spin  orbital  basis 

tl^  = <f>W 


(3.24) 

f W*  W'  \ 

(-/-•  = I,) 

It  should  be  noted  that  the  bullet  and  open  circle  notation  lose  their  meaning  of  occupied 
and  virtual  in  the  atomic  basis.  However,  we  are  still  using  occupied  and  virtual  for 
bullet  and  open  circle  even  in  the  atomic  basis  representation.  Using  the  fact  that  the 
only  transformation  that  will  have  any  effect  in  a determinantal  state  is  the  one  that  mixes 
unoccupied  orbitals  into  the  occupied  orbitals,  we  can  have  that  W"  = 0.  This  means 
that  the  reference  spin  orbitals 


= (f>*W'  + <I>°W° 


(3.25) 


are  constructed  by  first  orthonormalizing  the  bullet  (occupied)  atomic  orbitals  among 
themselves,  and  then  orthonormalizing  the  open  circle  (unoccupied)  atomic  orbitals  to 
the  bullet  (occupied)  space  and  among  themselves.  This  has  the  effect  that  the  occupied 
(bullet)  space  is  the  same  whether  defined  in  the  atomic  or  the  spin  orbital  basis.  That 
is  the  reason  why  we  keep  using  occupied  and  virtual  for  the  bullet  and  open  circle 
spaces  even  in  the  atomic  basis  representation.  We  have  then  created  a recipe  of  how 
to  construct  the  reference  state  from  the  atomic  orbitals  basis.  The  construction  of  the 
spin  orbitals  involves  what  is  called,  in  electronic  structure  theory,  the  MO  (molecular 
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orbital)  transformation.  This  is  a very  demanding  step  in  terms  of  computational  effort, 
since  it  scales  as  > K^,  and  the  atomic  integral  calculation  scales  as  w K^.  As  a result, 
despite  the  fact  that  the  equations  for  the  metric  and  dynamical  variables  get  a little  bit 
more  complicated,  it  is  of  great  computational  advantage  to  work  in  the  atomic  basis. 

The  spin  orbitals  that  are  delocalized  over  the  entire  molecular  system  are  also  called 

molecular  orbitals.  We  denote  the  set  of  parameters  of  the  coherent  state  with  respect  to 

the  reference  (spin  orbital,  or  molecular)  basis  as  and  with  respect  to  the  atomic 

basis.  The  following  transformations  exist 

2“*  = -I- 

(3.26) 

Z^nol  ^ 

and  also  the  basis  field  operators  transform  as 

b = wU 


(3.27) 


where  {a, a’}  are  the  field  operators  in  the  atomic  basis  representation. 
The  coherent  state  is  thus  equal  to 


N 


K 


t=l 


>=^  + 1 


N 


K 


= q;“*  n 


«=1 


j=N+l 


(3.28) 


at\at 


= 1^.-) 
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and  as  before,  we  work  with  the  unnormalized  coherent  state 

k“r=nK+  E “pji 

»=i  \ 

Note  that  because  the  anticommutation  relations  among  are  not  canonical^  it  is 

not  possible  to  write  the  coherent  state  as  an  exponential  with  the  atomic  parameters  and 
However,  it  is  still  possible  to  express  the  dynamical  orbitals 

K 

Xi^^i+  4>jZji,  (l<i<N)  (3.31) 

j=N+l 


^ I vac)  (3.29) 


N 

Xj  = - XZ  4-V’.^  (A"  + 1 < i < K) 

1=1 

in  the  same  simple  form 

K 

Xi  = ^i+  X]  (1  < * < N) 

j=N+l 


(3.32) 


(3.33) 


where 


N 

Xj  = 4>}~Y.  {N  + l<j<  K) 

1=1 

V = - (A*  -I-  A'z)"^ 

= - (a*  + A'^)  (a'  + 


(3.34) 


(3.35) 


1.  the  creation  and  annihilation  field  operators  in  the  spin  orbital  basis  are  canonical  since 
[6, 6^].^.  = 1,  however,  in  the  atomic  basis  [a,  = A the  anticommutation  is  not  canonical, 

since  the  atomic  overlap  matrix  A has  elements 

A.J  = {4>i\<l>j)  = J 4>i{x)4>j{x)dx 

and  {<?!>,  A'  being  non-orthononmal  makes  A different  of  the  unit  matrix 


(3.30) 
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with  the  overlap  being  given  by 


S = S{z'*,E',P\z,R,P)  = {z'\R',P'\z,R,P) 


= det(A*  + a'  z + z'^  a"  + z'^A°z) 


(3.36) 


We  now  define  the  following  matrices 


A*  A'  /• 
A° 


^ = A*  + z^A'^  + a'  z + z^A°z 


(3.37) 


and 


A°  = (t.  =A»  + ..A'  + AV  + rAV 


(3.38) 


which  naturally  occur  in  many  expressions.  Note  that  A*  and  A°  are  dependent  upon 
{z,z^}  and  respectively. 

The  metric  C can  now  be  expressed  in  the  atomic  orbital  basis  as 


C..  = 


dHnS 


- = -(a-^(a'  + 2^A°))^.^(a-^(a'  + 2+A°)) 


nk 


(3.39) 


c...= 


d'^lnS 


= ((aV  + A->)a»->  (vA'  + A->))^^{A->)^„. 


and 


Cz'z' 


dHnS 


The  mixed  derivatives  dz*V r involved  in  Cr  are 


d 


= (a%^  + A°)a°-^(u 


(3.40) 


^ + + <3.41) 


(3.42) 
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And  the  double  derivative  with  respect  to  the  nuclear  coordinates  is  given  by 


Crr  = V g^\n{z'  ,Fi\z,R%i^^^R>^R 


= Tr 


A— ^(/* 


(3.43) 


/• 


•—1/ 


/• 

z 


where 


P-'M) 

We  have  introduced  the  notation  of  a vertical  bar  to  indicate  that  the  derivative  of  the 
overlap  matrix  is  to  be  taken  with  respect  to  the  R dependence  of  one  side  only. 

Now  that  explicit  expressions  for  the  metric  have  been  given,  we  turn  to  the 
expression  for  the  energy  derivatives  appearing  in  the  right  hand  side  of  the  equation 
of  motion.  In  order  to  obtain  explicit  expressions  for  the  derivatives  of  the  energy  it 
is  appropriate  to  define  the  one-particle  density  matrix  (also  called  1-density  or  charge- 
bond-order  matrix)  as  follows 


r = (^J)A-'(/*  2'). 

which  was  obtain  from  the  kernel 

Tji{z*,z)  = {z\zy  {z\blbj\z) 
and  has  the  familiar  block  form 


(3.45) 


(3.46) 


(3.47) 
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The  energy  derivatives  then  become 
dE{z\z) 


+ /°)(/i  + Tr(Va6;„6r)J 

= ^ (a'^i;+  + A°)  ( u P )F  ^ ^ A*" 

where  we  recognize  the  Fock  operator 

F = h + Tr(yoj,;aj,r)^ 

= h-\-  Tr([(aa|66)  — (a6|fea)]r)^  . 

The  energy  derivatives  with  respect  to  the  nuclear  coordinates  are 


ki 


1=1 


(3.48) 


(3.49) 


F{z*,R',P\z,R,P)  I 

l^ZkZie^(^Rk-Ri)  / \ / \ 

= -5  E + Tr(v^.Ar)  - Tr(v^-,Arftr)  (3.50) 


+ 5Tr(Tr(v^^Ki;,»r)  j)  - Tr('Tr(v^^Ary,i.,tr)  J 
and  since  the  basis  set  is  independent  of  the  velocity  (or  momentum)  of  the  nuclei  the 
energy  derivatives  with  respect  to  the  momenta  become 


VpE{z\R',P\z,R,P) 


= fL 

R’=R,P'=P  Mk 


(3.51) 


We  now  have  explicit  equations  for  all  terms  in  the  equations  of  motion 
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and  we  can  try  to  simplify  the  main  equations  in  such  a way  that  they  become  efficient 
to  implement  in  a computer  program.  Starting  with  the  equation 


+ iCRR 


dE 

dz* 


(3.53) 


we  can  use  the  expressions  for  Cg>z,  Cr,  and  dEjdz*  to  further  reduce  it  to 


/• 

2: 


= ( -2  1°  )F 


/• 


(3.54) 


The  right  hand  side  has  a very  simple  interpretation:  the  Fock  operator  acts  on  the 
occupied  states  and  is  then  projected  on  the  virtual  space.  Thus  only  the  projection  on 
the  virtual  states  gives  rise  to  a change  in  the  2 coefficients.  This  is  also  the  way  that 
this  particular  equation  of  motion  is  implemented  since  it  does  not  involve  any  matrix 
inversion,  which  would  be  necessary  if  the  metric  did  not  have  this  particular  structure. 

In  the  equation  of  motion  that  involves  the  derivatives  of  the  energy  with  respect  to 
the  nuclear  coordinates, 


E{z\H',P\z,R,P)  I 


= 4E 


^ ZtZ,e‘‘{Rt  - Rl) 


/=! 
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+ Tr(v^-/r)-Tr(v^^Ar/.r)  (3.55) 
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we  can  see  that  it  will  be  necessary  to  compute  the  derivative  of  the  1-density  matrix 

r=  (^*j(^A*  + A'z  + z^A'^  + z^A°zy\r  z^) 


_M  1a— !(/. 


/• 
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V(/*  \i'  ^') 


with  respect  to  the  nuclear  coordinates,  which  can  obtained  as  follows 


= ; r 


)(V^.A-)(/' 


r\..- 
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(/•  .')(v^.a)(^Qa-‘(/'  ^') 


(3.56) 


(3.57) 


= r(v^,A)r. 

As  a result,  in  order  to  compute  the  energy  gradient  with  respect  to  the  nuclear  coordinates 
all  we  need  are  the  gradients  of  the  one-  and  two-electron  integrals  (which  are  available 
from  many  quantum  chemical  program  package)  in  addition  to  the  derivatives  of  the 
atomic  overlap  matrix. 

The  NDDO  Approximation 


The  NDDO  (neglect  of  diatomic  differential  overlap)  approach  employs  the  following 


approximations  [4]: 


32 


1.  Only  the  valence  electrons  are  considered  explicitly  in  the  calculations.  The  nucleus 
and  inner  shell  electrons  are  replaced  by  a fixed  core  function.  This  is  the  so  called 
core  approximation  [44]. 

2.  A minimum  basis  set  is  used.  This  basis  set  is  constructed  with  Slater  type  orbitals 

(STO’s).  The  advantages  of  using  STO’s  will  be  clarified  later.  The  STO’s  are 
products  of  a radial  function  and  a normalized  real  spherical  harmonics 

with  quantum  numbers  n,  /,  m, 

Xnlm  = Rnl{r)Yi^{e,(p), 


(or  ,'\n+l/2  (3.58) 

R ,(r)  = -n-1  -Cnir 

' [(2n)!]V2  ' • 

Here  Cn/  is  the  orbital  exponent  of  the  Slater  atomic  orbital. 

3.  The  overlap  integrals  are  neglected,  except  when  they  appear  in  the  one-electron 
resonance  integrals  (one-electron  two-center  integrals).  As  a result,  the  NDDO 
approximation  consists  of  applying  the  following  relation 


(1)  ^ SabX^{1)x^{1)  . (3.59) 


Since  we  are  treating  only  valence  electrons  and  using  STO’s  we  have  that 


(3.60) 


the  atomic  overlap  matrix  is  approximated  by  the  identity  matrix. 


4.  Two-electron  integrals  then  become 


^^C<Bc(#^(l)tf(2)|^^(l)^f(2))  (3.61) 
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We  have  chosen  the  NDDO  model  [45,  4]  since  it  is  considered  the  lowest  level 
at  which  there  exist  a basis  set  for  which  the  zero  differential  overlap  (ZDO)  [46] 
approximation  is  valid,  or  almost  valid  [47,  48].  The  NDDO  model  is  rotationally 
invariant  so  that  it  is  not  necessary  to  perform  any  spherical  averaging  of  the  atomic 
orbitals  with  angular  momentum  larger  than  zero.  As  a result,  unlike  the  other  ZDO 
models,  such  as  CNDO  (complete  neglect  of  differential  overlap)  and  INDO  (intermediate 
neglect  of  differential  overlap),  the  NDDO  model  includes  orbital  anisotropies.  In 
addition,  in  the  HF-NDDO  method  the  computational  dominant  step  is  the  diagonalization 
of  the  Fock  matrix. 

We  now  present  some  justifications  for  the  approximations  adopted  and  their  advan- 
tages: 

(1)  The  core  approximation  in  the  NDDO  model  is  reasonable  because  the  inner 
electrons  are  tighdy  bound  and  are,  therefore,  unlikely  to  be  significantly  perturbed  by 
changes  in  the  valence  shell.  Consequently,  most  chemical  processes  can  be  accurately 
described  by  valence  electrons  only.  Also,  there  are  solid  theoretical  foundations  for  the 
core  approximation  in  terms  of  an  effective  Hamiltonian  [49]. 

(2)  The  use  of  STO’s  as  building  blocks  for  the  molecular  orbitals  ensures  that  only 
a few  of  these  STO’s  are  necessary  to  give  a good  description  of  neutral  molecules 
or  positive  ions.  This  justifies  the  approximation  of  minimal  basis  set.  Only  in  the 
cases  of  anions  or  hypervalent  compounds  must  the  atomic  basis  set  be  more  flexible  to 
accommodate  the  extra  charge  density  or  extra  bond(s).  Adding  diffuse  or  polarization 
atomic  orbitals  is  one  way  of  solving  this  problem.  Another  is  to  allow  the  orbital 
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exponents  to  vary  with  the  atomic  charge.  This  provides  an  appropriate  description  for 
molecular  anions  [50],  and  also  allows  us  to  keep  the  minimum  basis  set  formalism.  The 
use  of  STO’s  is  useful  in  the  gradient  calculation.  That  is,  the  differentiation  of  an  STO 
function  with  respect  to  the  nuclear  coordinate  does  not  yield  higher  angular  momentum 
functions.  In  other  words,  once  the  basis  set  is  defined  in  terms  of  STO’s  no  higher 
angular  momentum  functions  are  needed  in  the  gradients  of  the  molecular  integrals. 

(3)  The  neglect  of  overlap  is  consistent  with  the  approximation  employed  in  the  two- 
electron  integrals.  However,  keeping  the  overlap  in  the  one-electron  resonance  integrals 
(one-electron  two-center  integrals)  is  important  since  these  integrals  are  responsible  for 
the  bonding  in  molecules.  That  is,  the  resonance  integrals 

= (x„(l)l(-iv?-^:?i)|x.(l)),  (3-62) 

which  involve  an  overlap  distribution  are  not  neglected. 

(4)  It  can  be  shown  [51,  52]  that  the  approximation  of  the  two-electron  integrals  by 
the  neglect  of  diatomic  differential  orbital  is  correct  up  to  second  order  in  the  overlap 
expansion.  It  is  also  possible  to  evaluate  the  three-  and  four-center  integrals  in  terms 
of  overlap  and  two-center  integrals  [52].  However,  if  we  parametrize  the  two-electron 
integrals  with  respect  to  experimental  data  the  errors  caused  by  neglecting  the  three-  and 
four-center  integrals  can  be  minimized. 

The  END-SD(NDDO)-FGWP  Method 


Since  the  dynamical  equations  are  dependent  upon  the  norm  we  start  by  applying 
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the  NDDO  approximation  to 


S = det(A*  + A' z + z'^A"  -f  z'^A^z). 
The  NDDO  approximation  implies  that 


A*  = r,  A'  = A"  = o,  A°  = r 

so  that  the  overlap  takes  the  same  form  as  in  the  orthonormal  basis, 

S = det(7*  + z'^z) 


as  do  the  A matrices 

A*  = /•  + zh 


A°  = 7°  + zzK 


Also,  the  density  matrix  is  simplified  to 


r=  1^*1  A— ^(7*  z^)  = 


(^'’yr  + z'z)-'{r  zt) 


(3.63) 


(3.64) 


(3.65) 


(3.66) 


(3.67) 


and 

u = -(A'^  + A°2)(A*  + A'z)-^  = -z 


= -(A*  + z^A'^)-^(A'  + z^A®)  = -z+ 

As  a result,  the  relevant  terms  for  the  dynamical  equations  become 

Metric 


d'^lTiS 


-((r  + z*z)-'z^).,{{r  + z<z)-'z\, 


(3.68) 


(3.69) 
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(3.70) 


Im  ^ki 


(3.71) 


d 


= (/°  + ^zt)-l(_^  /'>)V^-^A|i:  (/*  + zM 


K\-i 


Crr  = V ^\n{z'  ,Ii\z,R)\^'^^^Il-^R 


(3.72) 


= Tr 


(/•  + 2^2)-^(/* 


V^AI  ^ 


-(r  + z'z)-'(r  z^)v^A\(''){r  + z'z)-\r  z«)v^a/^‘ 


Energy  Gradients 


dE{z\z) 


8z* 


= ((/°  + 22')-‘(-z  r)F{’'\r  + zU)-' 


ki 


V^E{z\B!,P\z,R,P)\^ 

"*  \R'=R,P>= 


^ZkZ,e^(^Rk-Ri) 


1 

=—y 

2 ^ 


1=1 

l^k 


- RiP 


+ Tr(Vjj,Ar)-Tr( 


V^AThr) 


+ 5Tr('Tr(v^^Ki..ir)  - Tr(Tr(v^^ArK.,,.ir)  F 


(3.73) 


(3.74) 


(3.75) 


Vpy{z\R',P\z,R,P) 


= A 

R'=R,P'=P  Mk 


(3.76) 
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where  the  Fock  operator  must  also  be  considered  under  the  NDDO  approximation 

f = A + Tr(v;^S'>°r)^ 

(3.77) 

= h + TT{[{aAaA\bBbB)  ~ {aAbA\bBaB)]^)a  • 

Caution  must  be  exercised  when  considering  the  derivatives  of  the  atomic  orbital 
overlap  A with  respect  to  the  nuclear  position.  If  we  make  the  NDDO  approximation 
for  the  overlap,  that  is,  A = /,  then  take  the  derivative  we  obtain  no  contribution. 
As  a result,  both  Cr  and  Crr  vanish  and  we  do  not  have  any  coupling  between  the 
electrons  and  the  nuclei.  This  is  not  consistent  with  the  ab  initio  END  equations,  since  the 
NDDO  approximation  should  not  affect  the  proper  coupling  between  electrons  and  nuclei. 
Also,  making  the  NDDO  approximation  before  taking  the  derivatives  is  inconsistent  with 
the  analytical  gradient  method.  That  is,  the  contribution  of  the  two-center  one-electron 
(resonance)  integrals  to  the  energy  gradient  would  then  vanish,  which  is  not  true  when 
finite  difference  method  is  employed.  Consequently,  we  should  first  take  the  derivative  of 
the  non-approximated  atomic  orbital  overlap  and  then  perform  the  NDDO  approximation 
on  the  resulting  expression. 

AMI  Model 

The  AMI  (Austin  Model  1)  approach  to  the  NDDO  approximation  [53]  has  been 
extensively  tested  and  has  been  shown  to  reproduce  geometry,  heats  of  formation,  dipole 
moments,  and  ionization  potentials  very  well  [53-55].  In  addition,  the  overall  perfor- 
mance of  the  AMI  for  computing  molecular  properties  (excitation  energies,  vibrational 
frequencies  and  intensities)  and  quantities  relevant  for  chemical  reactions  (enthalpy  of 
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reaction  and  activation  energies)  is  also  satisfactory  [55].  It,  thus,  seems  appropriate  to 
use  the  AMI  realizations  of  the  NDDO  method  in  our  END-SD(NDDO)-FGWP  model. 
Results  for  dynamical  properties  like  intra-  and  intermolecular  charge  transfer,  transition 
state  spectroscopy,  photodissociation,  etc,  can  be  used  to  check  the  performance  of  the 
END- AMI  method.  In  case  this  approach  fails  a different  realization  of  the  NDDO 
method  can  implemented  where  the  parametrization  should  also  include  experimental 
dynamical  properties  in  the  data  base. 

Before  we  present  the  AMI  model  we  need  to  consider  first  its  Fock  matrix  elements 
[56,  53,  55]: 

1 . Diagonal  terms 


R 


= u, 


(tA/iA 


Bi^A 


(3.78) 


B^A  Xb 


2.  Off-diagonal  terms  on  the  same  atom 


Bi^A 


+ l:P^^Al'A[^y'Al'A\^^AVA)  - {tiAiiA\vAt'A)\ 


(3.79) 


+ '^'Y^P\B<TBit^AVA\^B(^B)- 


B^A  Xb  ^b 
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3.  General  terms  on  different  atoms 


I'A  <TB 


where  <^b)  are  atomic  orbitals  centered  in  A {B)  and  P^x  is  the  bond 

order  matrix. 

In  order  to  complete  the  definition  of  the  AMI  model  we  need  expressions  or  numbers 
for  each  one  of  the  terms  that  go  into  the  elements  of  the  Fock  matrix.  That  is, 

!•  UfiAHA  one-electron  one-center  terms  and  are  parametrized  from  spectroscopic 
data  for  valence  states  of  atom  A and  its  ions; 

2.  Za  is  the  atomic  number  of  atom  A minus  the  number  of  core  electrons; 

3.  The  nuclear  repulsion  energy  is  given  by 


where  af,  bf,  and  cf  are  parameters; 

4.  The  resonance  integral  is  made  proportional  to  the  atomic  orbital  overlap.  The 
proportionality  constant  is  the  average  of  the  parameters; 


Eab  = ZaZb{sasa\sbsb) 


4 


(3.81) 


5.  The  repulsion  integrals  are  expanded  in  terms  of  multipole  Mi^  and  are  approximated 
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by  the  following  expression 

00  OO  imtn 

/l=0  /2— 0 


where, 


Rij  = ^AB  + + vf,D^  (3.83) 

with  tj  = ...,±2,±1,0  depending  upon  the  order  of  the  dipole,  p and  D are  derived 
parameters,  that  is,  D = /^(C/i)  and  p = p{U,i„,D).  In  the  case  of  a minimum  sp 
basis  set,  there  are  22  distinct  repulsion  integrals.  The  inclusion  of  d orbitals  increase 
this  number  to  450,  which  make  this  type  of  approximation  for  the  repulsion  integrals 
cumbersome  to  be  extended  beyond  sp  basis. 


The  number  of  parameters  in  the  AMI  model  ranges  from  13  to  16.  Around  150 
molecules  were  included  in  the  reference  data  base  used  to  obtain  those  parameters.  These 
numbers  can  give  an  idea  of  how  difficult  is  to  obtain  parameters  for  a semiempirical 
method. 

Implementation 

The  END-SD-FGWP  method  has  been  implemented  in  the  ENDyne  program  package. 
This  program  is  capable  of  performing  simultaneous  optimization  of  the  electronic  and 
nuclear  (geometry,  for  classical  nuclei)  wave  functions.  ENDyne  also  can  perform  the 
time  evolution  of  a molecular  system  state.  The  input  consists  of  the  initial  states 
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(electronic  and  nuclear),  number  and  type  of  particles  (electrons,  nuclei,  nuclear  masses, 
atomic  number),  initial  and  final  time  of  the  evolution,  the  (Gaussian)  basis  set  for  each 
nucleus,  and  the  type  of  differential  equation  solver.  The  program  works  in  the  laboratory 
frame  and  uses  atomic  units.  As  output,  ENDyne  generates  a file,  called  restart  which 
contains  the  history  (the  2 parameters  and  phase)  of  the  system  evolving  in  time.  Another 
program,  called  EVOLVE,  is  used  to  extract,  plot,  and  manipulate  relevant  information 
like  nuclear  position,  electronic  and  nuclear  momenta,  phase,  electron  density  and  atomic 
population,  etc,  from  the  restart  file. 

A functional  diagram  of  the  ENDyne  program  can  be  seen  in  the  Figure  3-1. 
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Figure  3-1:  Functional  diagram  of  the  ENDyne  program. 


The  main  program  endyne  calls  getinp  and  decides  the  structure  and  type  of  calcula- 
tion. After  that,  it  makes  a dummy  run  and  determines  exactly  how  much  memory  will  be 
needed  for  the  calculation.  Then,  the  rundyn  is  called  and  it  takes  control  of  the  ENDyne 
calculation.  It  calls  dynini  to  initialize  the  calculation  for  the  models  (quartic,  Morse, 
and  oscillator  potentials)  and  theories  (Hartree-Fock,  AGP,  RPA).  It  then  decides  if  the 
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calculation  is  an  optimization  (frprnm,  djpmin),  optimization  (Jrprmn,  dfpmin)  followed 
by  evolution  {deint,  drint,  odeint),  only  evolution(dc,  driveb,  odeint),  or  exact  oscillator 
(xctosc).  Dynfun  is  then  called  and  computes  the  gradient  of  the  energy  with  respect  to 
the  wave  function  parameters  in  the  case  of  optimization  or  the  derivatives  of  all  wave 
function  parameters  with  respect  to  time  for  an  evolution  calculation.  The  appropriate 
forces  (gradients)  for  the  semiempirical  or  ab  initio  model  is  then  called.  In  the  case 
of  the  Hartree-Fock  theory  (END-SG-FGWP)  the  atomic  integrals  and  overlap  and  their 
derivatives  are  computed  every  time  the  force  routine  {hffrc)  is  called.  The  onedrv  and 
twoint  routines  are  drivers  for  one-  and  two-electron  integral  calculations,  with  averll, 
averll,  and  disint  being  the  interfaces  with  the  ABACUS  integral  package. 

In  the  case  of  AMI  model  the  structure  is  the  same,  except  that  the  drivers  onedrv 
and  twoint  calls  some  other  appropriate  interfaces,  instead  of  calling  averll,  averll, 
and  disint.  In  addition,  the  implementation  of  the  AMI  method,  necessitates  the 
development  of  equations  and  source  code  for  computing  analytical  derivatives  of  the 
integrals  (repulsion  and  overlap)  and  of  the  energy.  The  molecular  orbital  programs 
AMPAC  2.1  and  MOPAC  6.0  based  upon  the  AMI  model  do  not  perform  analytical 
derivatives  at  the  SCF  (self-consistent  field)  level.  However,  analytical  derivatives  for 
similar  expressions  of  the  repulsion  integrals  have  already  been  performed  [57],  which 
insures  that  the  AMI  model  implementation  is  viable. 


CHAPTER  4 

H+  + H,  He,  AND  H2  CHARGE  TRANSFER  COLLISION 


In  this  Chapter  we  present  results  for  the  electron-nuclear  dynamics  that  takes  place 
when  H"^  collides  with  H,  He,  or  H2. 

Such  a study  consists  of  computing  the  probability  of  the  charge  transfer,  energy 
transfer,  etc  as  a function  of  the  impact  parameter  and,  in  the  case  of  the  H2  target,  also 
molecular  orientation.  Plots  of  the  transfer  probability,  vibrational  excitation,  and  cross 
sections  versus  the  scattering  angle  are  presented  and  comparisons  with  other  theoretical 
approaches  as  well  as  experimental  data  discussed. 

This  chapter  has  the  following  subsections:  the  H'*'  + H,  + He,  and  + H2  results,  with 
a fourth  subsection  containing  general  conclusions. 

Before  presenting  the  results  some  general  considerations  are  in  order. 

i)  All  the  calculations  were  performed  with  the  ENDyne  program  package  using  the 
following  computers:  IBM  RS/6000-550,  Sun  4/490,  Sun  690MP,  and  Cray  Y-MP/432. 

ii)  An  Is2s2p  basis  set  was  used  for  both  H and  He.  This  basis  sets,  namely  pVDZ 
for  H:  [4slp]/(2slp)  [58]  and  6-31G*  for  He:  [4slp]/(2slp)  [59]  are  listed  in  Table  4-1. 
In  the  pVDZ  basis  set  the  s-exponents  were  scaled  by  1.2^  = 1.44,  as  was  done  in  the 
original  DZP  basis  set  by  Dunning  [60]. 

iii)  Some  of  physical  properties  for  the  atoms  and  molecules  involved  in  this  work  are 
presented  in  Table  4-2. 

iv)  Due  to  the  physics  of  the  problem  (a  proton  hitting  atomic  and  molecular  targets)  the 
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initial  wave  function  consisted  of  a totally  localized  electron  density  on  the  H,  He,  and 
H2  targets.  In  other  word,  the  starting  wave  functions  are  not  the  ground  states  of  the 
complex  species  (H-H)"^,  (H-He)"^,  or  (H-H2)‘^.  It  should  be  mentioned,  however,  that 
the  atomic  and  molecular  targets  are  in  their  electronic  (and  vibrational)  ground  states. 
The  difference  in  energy  associated  with  these  initial  wave  functions  and  the  ground 
state  energy  of  the  (super)molecule  ion  species  is  due  to  the  localization  of  the  electronic 
density  around  the  targets. 


Table  4-1:  Contraction  coefficients  c and  exponents  a for  the  basis  sets  used  in  this  work. 


H/pVDZ“ 

He/6-3 

c 

a 

c 

a 

Is  orbital 

0.01969 

18.7344 

0.02311 

38.4216 

0.13798 

2.82528 

0.15468 

5.77803 

0.47815 

0.64022 

0.46930 

1.24177 

2s  orbital 

0.50124 

0.17568 

0.29796 

1.00000 

2p  orbital 

1.00000 

0.72700 

1.10000 

1.00000 

a:  Ref.  [58],  b:  Ref.  [59]. 


Table  4-2:  Ionization  potentials,  electron  affinities  and  polarizabilities  of  H,  He  and  H2. 


Ionization  Potential 
(eV) 

Electron  Affitinity 
(eV) 

Polarizability 
(lO'^"*  cm^) 

H* 

13.598 

0.754 

0.6668 

He® 

24.587 

not  stable 

0.2050 

H2® 

15.427 

— 

0.803 

a:  Ref.  [61]. 
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The  + H Collision 

The  study  of  H'*'  + H collisions  is  of  special  interest  because  of  its  great  simplicity 
but  rich  physics.  This  has  inspired  much  theoretical  and  experimental  work  involving 
charge  transfer  (total  and  differential  cross  sections),  excitation,  alignment  measurements, 
etc.  As  a result,  H'*’  + H collisions  provide  a good  test  for  any  time-dependent  method 
that  intends  to  describe  electron  and  nuclear  dynamics.  It  should  be  noted  that  since  this 
is  a one  electron  system  the  effects  of  electronic  correlation  are  absent.  Consequently, 
the  same  methods  should  be  tested  for  systems  with  more  than  one  electron,  since  like 
in  stationary  quantum  chemistry,  the  electronic  correlation  effects  can  be  essential  when 
describing  molecular  properties. 

Most  of  the  theoretical  calculations  on  the  H'*’  + H system  are  done  by  close- 
coupling methods  [25,  24].  Since  these  methods  use  prescribed  trajectories  they  are 
most  commonly  applied  to  collision  energies  above  1 keV  in  the  laboratory  frame.  One 
of  the  investigations  in  this  chapter  includes  the  comparison  of  dynamical  trajectories 
obtained  with  the  END-SD-FGWP  method  with  prescribed  Coulomb  trajectories  [34]. 
Below  1 keV  there  are  just  a few  theoretical  studies  where  the  PSS  (perturbed  stationary 
states)  method  is  used  to  separate  the  nuclear  and  electronic  degrees  of  freedom  [62]. 

Total  Cross  Sections 

The  scattering  of  a proton  by  a hydrogen  atom  has  already  been  studied  by  the 
END-SD-FGWP  method  [38].  In  that  work  an  Is2s2p  basis  set  was  also  employed  but 
each  function  was  expressed  in  terms  of  six  primitive  Gaussian  functions.  Since  in  this 
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work  we  are  using  a smaller  number  of  primitives  we  decided  to  compute  the  total  cross 
sections  at  several  energies  and  compare  them  with  experimental  and  previous  END 
work.  This  comparison  should  provide  us  with  some  idea  of  how  the  basis  set  affects 
the  total  cross  section. 

We  have  used  the  following  initial  conditions:  i)  the  hydrogen  atom  is  in  its  ground 
state  with  the  following  coordinates  (20.0,  0.0,  0.0)  in  a.u.  with  respect  to  the  laboratory 
frame;  ii)  the  proton  is  placed  at  the  (0.0,  b,  0.0)  coordinates,  with  b (the  impact  parameter) 
varying  from  0.0  to  15.0  a.u.;  iii)  the  hydrogen  atom  has  zero  momentum;  iv)  the  proton 
has  initial  momentum  corresponding  to  the  energy  of  the  collision;  v)  the  time  evolution 
of  the  system  is  carried  out  until  the  nuclei  are  separated  by  more  than  1(X)  a.u.  The 
transition  probability  is  obtained  by  projecting  the  asymptotic  state  of  the  proton  on 
the  eigenstates  of  the  moving  proton  for  elastic  processes  or  hydrogen  atom  for  charge 
exchange  [38]. 

Total  cross  sections  of  the  scattering  of  H'*’  on  H are  obtained  by  numerically 
integrating  [63,  64] 

CX) 

<7{E)  = 27t  / bP{b;E)db  (4.1) 

0 

for  several  energies  E (0.5,  1.5,  5.0  keV,  e.g.),  where  b is  the  impact  parameter  and 
P{b]  E)  the  transfer  probability  as  a function  of  P(6;  E)  for  a given  E.  The  upper 
integration  limit  is  chosen  to  be  between  11.0  and  12.0  a.u..  This  choice  is  reasonable, 
since  the  transfer  probability  in  this  interval  is  smaller  than  10~^. 

For  a collision  energy  of  0.5  keV,  we  obtain  the  results  summarized  in  the  Table 
4-3,  where  also  the  total  cross  section  is  shown  as  calculated  using  different  integration 
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steps  (increments  in  the  impact  parameter).  As  we  can  notice  from  Table  4-3,  the  best 
compromise  between  intervals  in  the  impact  parameter  and  the  cross  section  accuracy 
is  the  last  row,  where  we  have  an  increment  of  0.10  a.u.  in  the  oscillatory  region  (0.0 
- 6.0  a.u.)  and  0.50  a.u.  in  the  monotonically  decreasing  region  (6.0  - 11.0  a.u.),  see 
Figure  4-2  for  visualization. 


Table  4-3:  Total  cross  section  for  + H collision  at  500  eV  using  pVDZ  basis  set. 


Interval  (a.u.) 

Increment  in  b 
(a.u.) 

Interval  (a.u.) 

Increment  in  b 
(a.u.) 

Cross  Section 
([a.u.]2) 

0.0  - 11.0 

0.05 

64.0583 

O 

• 

1 

O 

• 

o 

0.10 

64.0583 

0.0  - 11.0 

0.15 

64.0671 

0.0  - 6.0 

0.05 

6.0  - 11.0 

0.15 

64.0589 

0.0  - 6.0 

0.10 

6.0  - 11.0 

0.50 

64.0521 
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Impact  Parameter  (a.u.) 

Figure  4-2:  Weighted  transition  probabilities  for  total  electron  transfer  at  500  eV  as  a function 
of  impact  parameter  from  ENDyne.  All  data  in  atomic  units. 

Table  4-4  summarizes  the  total  cross  section  results  obtained  by  the  END-SD-FGWP 
method  using  two  different  sets  of  primitive  functions  in  the  basis  sets. 


Table  4-4:  Total  transfer  cross  sections  for  H*  colliding  with  a H (xlO"*^  cm^). 


Total  Transfer  Cross-section. 

Collision 

END-SD-FGWP  Theory 

Experiment^*’ 

Energy  (eV) 

pVDZ  basis 

STO-6G  basis 

10 

35.93 

36.37 

37.0+5.7 

100 

24.31 

25.60 

23.7+3.5 

500 

17.94 

19.44 

18.9+3.2 

1000 

16.55 

16.78 

16.3±2.9 

a:  Newman  et  ai.  Ref.  [65],  b:  Gealy  and  Van  Zyl,  Ref.  [66]. 


50 


As  we  can  see  from  Table  4-4  the  basis  set  do  not  have  a large  effect  in  the  calculated 
integral  cross  section  for  electron  transfer  in  the  + H collision.  In  addition,  as  noted 
in  early  studies  the  END-SD-FGWP  method  provides  excellent  results  for  charge  transfer 
cross  sections  when  compared  with  experimental  data. 

Differential  Cross  Sections 

Since  the  integral  cross  section  involves  an  integration  much  of  the  detailed  behavior 
of  a charge  exchange  collision  is  masked.  Such  details  can  be  obtained  from  state-to- 
state  cross  section,  differential  cross  section,  and  transfer  probability  as  a function  of  the 
scattering  angle  and  consequently  should  be  the  basis  for  theoretical  model  analysis.  We 
present  results  for  the  transfer  probability  and  the  reduced  differential  cross  section  as 
a function  of  the  scattering  angle  for  the  H'*'  colliding  with  H at  several  energies  (0.25, 
0.41,  0.50,  0.70,  1.0  keV). 


The  expression  for  the  differential  cross  section  (aj)  at  a given  collision  energy  E 


is  [63,  64] 


CTd  = 


da{e)  b P{b) 


dO, 


sin  9 


d6 

IE 


(4.2) 


where  b,  6,  P{b),  dil  = sin  0 dO  d(f>  are  the  impact  parameter,  the  scattering  angle,  the 
transfer  probability,  and  the  solid  angle,  respectively.  The  modulus  of  db/d6  is  used  to 
insure  that  the  differential  cross  section  is  always  a positive  quantity.  In  addition,  since 
the  scattering  angle  by  definition  always  lies  between  zero  and  tt,  we  have  that  sin  0 > 0. 
It  should  be  noted  that  the  differential  cross  section  depends  upon  the  derivative  and  is 


not  integrated  over  any  variable  so  that  more  details  of  the  collision  process  are  revealed. 
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We  should  now  make  the  connection  between  the  classical  deflection  function  0 and 
the  scattering  angle  d.  The  classical  deflection  function  0 is  given  by  [63,  64] 


where  V(r)  is  the  central  potential,  E the  collision  energy,  and  the  distance  of  closest 
approach  a is  obtained  by  solving  the  equation 


(4.4) 


From  the  expression  for  the  deflection  function  we  can  see  that  0 may  have  any  value 
between  — oo  and  x,  and  as  we  mentioned  before  0 < 0 < tt.  As  a result,  the  two  angles 
have  to  be  connected  by 


Q = TjO  — 2nv 


where  t;  = ±1  and  n is  a positive  integer  (or  zero).  For  a given  0,  the  numbers  t)  and 
n are  chosen  so  that  6 lies  between  zero  and  tt.  We  can  also  show  that 


de  dd 

de 

de 

db  ~"^db  ^ 

db 

— 

db 

(4.6) 


and  consequently. 


(4.7) 
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Trajectory  Analysis 


As  mentioned  before  time-dependent  methods  using  prescribed  trajectories  are  not 
appropriate  for  calculating  the  differential  cross  sections  since  they  involve  the  derivative 
of  the  deflection  function  0 with  respect  to  the  impact  parameter  b.  The  END-SD-FGWP 
method  with  its  dynamical  trajectories  should  be  more  suitable.  We  present  a comparison 
of  0 versus  6 using  bare  nuclei  and  a hydrogen  atom  potential.  The  central  potential  in 
the  case  of  bare  nuclei  has  the  following  form 

r)  = z=  r~^  (4.8) 

r 


for  two  protons.  The  resulting  analytical  solution  for  the  relation  between  impact 
parameter  and  scattering  angle  is  the  so-called  Rutherford  scattering  formula 


This  is  certainly  not  a realistic  interaction  potential  between  a proton  and  a hydrogen 
atom.  Instead,  using  the  exact  wave  function  for  the  ground  state  of  a hydrogen  atom 
we  obtain  the  following  potential  for  the  electrostatic  interaction  between  a proton  and 
a hydrogen  atom  in  the  ground  state  [67] 


I/(r)  = (l-|-r-^)e-2". 


(4.10) 


There  is  no  analytical  solution  of  the  classical  deflection  function  for  this  potential,  so 
we  have  to  solve  it  numerically  and  the  results  are  plotted  in  Figure  4-3. 
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Figure  4-3:  Scattering  angle  as  a function  of  impact  parameter  using  ENDyne,  bare  nuclei,  and 
hydrogen  atom  potential  for  the  + H system.  Eneigy:  500  eV.  Basis  set:  pVDZ.  Scattering 
angle  in  degrees  and  impact  parameter  in  a.u..  Full  line:  ENDyne,  dotted  line:  bare  nuclei,  and 
dashed  line:  hydrogen  atom  potential. 

A more  detailed  illustration  of  the  differences  between  these  three  trajectories  can 
be  seen  in  Figure  4-4.  It  is  clear  from  Figure  4-4  that  the  errors  due  to  prescribed 
trajectories  can  be  quite  large  for  small  impact  parameter.  We  also  notice  that  the 
dynamical  trajectory  lies  between  the  bare  nuclear  and  the  hydrogen  atom  potentials. 
This  is  due  to  the  fact  that  in  the  H'*’  + H system  charge  transfer  is  expected  so  the 
interaction  potential  should  be  represented  by  a combination  of  the  bare  nuclei  and  the 
atomic  potential.  The  contribution  of  each  extreme  potential  (bare  nuclei  and  atomic) 
is  related  to  the  amount  of  charge  transfer  during  the  collision.  Consequently,  the  use 
of  prescribed  trajectories  in  calculating  properties  that  depends  upon  the  scattering  angle 
becomes  doubtful  for  energies  below  1 keV  and  for  small  impact  parameters. 
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Figure  4-4:  Scattering  angle  as  a function  of  impact  parameter  using  ENDyne,  bare  nuclei,  and 
hydrogen  atom  potential  for  + H system.  Small  impact  parameter  region.  Energy:  500  eV. 
Basis  set:  pVDZ.  Scattering  angle  in  degrees  and  impact  parameter  in  a.u..  Full  line:  ENDyne, 
dotted  line:  bare  nuclei,  and  dashed  line:  hydrogen  atom  potential. 

Transfer  Probability 


As  we  can  see  from  the  expression  for  the  differential  cross  section,  equation  (4.7), 
the  transfer  probability  is  an  important  contributing  factor  for  this  quantity.  We  then 
present  in  this  subsection  END-SD-FGWP  results  for  this  probability  as  a funciton  of  the 
scattering  angle  and  compare  them  with  the  experimental  data. 

Before  performing  these  comparisons  we  present  a brief  discussion  of  the  experi- 
mental procedure  and  try  to  identify  the  main  sources  of  experimental  uncertainties.  The 
experimental  procedure  consists  of  [68,  69]  a mass  analyzed  H"''  beam  issued  from  a 
single  discharge  ion  source  crossed  with  a thermal  H beam.  Ions  and  atoms  scattered 
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at  a given  angle  are  selected  and  analyzed.  The  overall  energy  resolution  varies  from 
0.8  eV  up  to  2 eV  in  the  250  eV  — 2000  eV  energy  range.  The  uncertainties  affecting 
the  experimental  results  arise  from  i)  the  finite  angular  resolution  AO  (±  0.07°  for  small 
angles  up  to  ± 0.2°  for  large  angles),  ii)  the  errors  in  the  dissociation  fraction  determina- 
tion (H/H2  ratio  in  the  collision  region),  and  iii)  the  accuracy  of  the  relative  calibrations 
of  the  cross  sections  for  the  different  processes  (elastic,  charge  transfer,  and  excitation). 
Because  the  rapid  variation  of  the  transfer  probability  and  the  differential  cross  section 
with  the  angle,  more  particles  arise  from  the  region  6 — AO  than  from  the  region  0 -f-  A^. 
Thus  the  center  angle  is  slightly  larger  than  the  one  corresponding  to  the  most  prob- 
able scattering.  However,  this  error  never  exceeds  0.1°  [68,  69].  As  we  shall  see  in 
the  comparison  between  the  experimental  data  and  the  END-SD-FGWP  results  the  finite 
angular  resolution  can  account  for  most  of  the  damping  in  the  oscillatory  behavior  of  the 
observed  versus  calculated  transfer  probability  and  differential  cross  sections.  Problems 
with  the  measured  direction  of  the  incident  beam  with  respect  to  the  detection  system 
can  be  eliminated  by  studying  the  scattering  on  both  sides  of  the  incident  beam  direction 
and  choosing  the  zero  index  from  the  symmetry  of  the  data  curves.  In  addition  to  these 
uncertainties,  the  experimental  results  are  normalized  to  the  theoretical  results  [70]  for 
the  elastic  differential  cross  section  at  small  angles. 

From  the  discussion  of  the  experimental  uncertainties  it  seems  clear  that  the  measured 
transfer  probability  versus  the  scattering  angle  is  less  susceptible  to  errors  than  the 
differential  cross  section  and  they  are  independent  of  scaling  factor.  Thus  in  Figures 
4-5-4-9  we  present  the  results  for  the  transfer  probability  calculated  by  the  ENDyne 
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program  and  compare  with  the  experimental  data.  The  theoretical  procedure  is  the 
same  as  described  in  the  calculation  of  the  integral  cross  section.  The  only  additional 
information  needed  is  the  asymptotic  momenta  of  the  proton,  which  are  used  to  compute 
the  scattering  angle,  that  is 

Px  Px  1 / .1  1 1 \ 

COS0  = = . , = = (4.11) 

Ptot  y/Px+Py  y/1  -f  {PylpxY 

where  px  and  py  are  the  x-  and  y-  components  of  the  asymptotic  proton  momentum. 


Scattering  Angle  (deg) 

Figure  4-5:  Transfer  probability  versus  the  scattering  angle  (in  degrees).  Comparison  between 
the  experimental  and  ENDyne  results  for  the  H*  + H system.  Energy:  250  eV.  Basis  set:  pVDZ. 
Experimental  angtilar  resolution  ±0.6°,  ’+’  ENDyne,  from  Ref.  [68],  and  ’x’  from  Ref.  [69]. 
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Figure  4-6:  Transfer  probability  versus  the  scattering  angle  (in  degrees).  Comparison  between 
the  experimental  and  ENDyne  results  for  the  + H system.  Energy:  410  eV.  Basis  set:  pVDZ. 
Experimental  angular  resolution  ±0.2°  — ±0.6°,  '+'  ENDyne,  firom  Ref.  [68],  and  ’x’  from 
Ref.  [69]. 


Figure  4-7:  Transfer  probability  versus  the  scattering  angle  (in  degrees)  for  the  H*  + H system. 
Energy:  500  eV.  Basis  set:  pVDZ.  ’+’  ENDyne. 
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Figure  4-8:  Transfer  probability  versus  the  scattering  angle  (in  degrees).  Comparison  between 
the  experimental  and  ENDyne  results  for  the  + H system.  Energy:  700  eV.  Basis  set:  pVDZ. 
Experimental  angular  resolution  ±0.02°  — ±0.6°,  ’+’  ENDyne  and  from  Ref.  [68]. 


Figure  4-9:  Transfer  probability  versus  the  scattering  angle  (in  degrees).  Comparison  between 
the  experimental  and  ENDyne  results  for  the  H"^  + H system.  Energy:  1000  eV.  Basis  set:  pVDZ. 
Experimental  angular  resolution  ±0.07°  — ±0.6°,  ’+’  ENDyne  and  from  Ref.  [68]. 
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Both  the  qualitative  and  quantitative  features  (position  of  maxima  and  minima)  of 
the  experimental  data  are  very  well  reproduced  by  the  END-SD-FGWP  method  The 
quantitative  agreement  can  be  made  much  better  if  we  apply  an  angular  resolution  window 
to  the  theoretical  data.  This  would  simulate  the  damping  in  the  experimental  data  and 
explain  why  the  oscillations  in  the  transfer  probability  do  not  extend  form  zero  to  unity. 
There  is  still  some  discrepancies,  mainly  at  high  energies  and  large  angles,  between  our 
theoretical  treatment  and  the  experimental  data.  We  believe  that  the  extra  features  in 
the  theoretical  results  are  real  and  were  not  picked  up  by  the  experiment  due  to  the  poor 
angular  resolution.  It  also  should  be  mentioned  that  the  semiclassical  treatment  [70]  using 
a three-term  molecular  basis  provides  results  in  perfect  agreement  with  the  experimental 
data.  This  is  puzzling  since  the  experimental  data  have  a damping  of  the  oscillation 
due  to  the  finite  angular  resolution  and  the  theoretical  work  does  not  mention  how  this 
damping  is  handled.  Because  of  this  the  quantitative  agreement  between  the  semiclassical 
method  and  the  experimental  results  seems  to  be  fortuitous  or  the  true  meaning  of  the 
theoretical  data  is  difficult  to  understand.  We  should  add  that  this  damping  is  also  not 
observed  in  other  theoretical  treatments  [34]. 

Reduced  Differential  Cross  Section 


The  experimental  results  available  for  the  H'*’  + H collision  are  often  presented  as 
’reduced’  differential  cross  sections  (p)  which  are  expressed  as  [68] 


P 


be  P{b) 


db 


(4.12) 
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The  reason  for  this  is  due  to  the  rapid  variation  of  the  differential  cross  sections  with 
the  scattering  angle.  So,  in  the  following  Figures  (4-10-4-14)  we  present  both  elastic 
and  charge  transfer  reduced  cross  sections  calculated  using  the  ENDyne  program  and 
compare  them  with  the  experimental  data. 


0.0  1.0  2.0  3.0  4.0  5.0  6.0  7.0  8.0 


Scattering  Angle,  0 (deg) 

Figure  4-10:  Reduced  differential  cross  sections  versus  the  scattering  angle.  Experimental  and 
the  ENDyne  results  for  the  H"^  -h  H system.  Eneigy:  250  eV.  Basis  set:  pVDZ,  ENDyne 
transfer,  ENDyne  elastic,  ’x’  transfer  and  ’o’  elastic  from  Ref.  [68]. 
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Figure  4-11:  Reduced  differential  cross  sections  versus  the  scattering  angle.  Experimental  and 
the  ENDyne  results  for  the  -i-  H system.  Energy:  410  eV.  Basis  set:  pVDZ,  ENDyne 
transfer,  ’+’  ENDyne  elastic,  ’x’  transfer  and  ’o’  elastic  from  Ref.  [68]. 


Scattering  Angle,  0 (deg) 

Figure  4-12:  Reduced  differential  cross  sections  versus  the  scattering  angle.  Experimental  and 
the  ENDyne  results  for  the  H system.  Energy:  500  eV.  Basis  set:  pVDZ,  ’*’  ENDyne 
transfer,  ’+’  ENDyne  elastic,  ’x’  transfer  and  ’o’  elastic  from  Ref.  [68]. 
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Figure  4-13:  Reduced  differential  cross  sections  versus  the  scattering  angle.  Experimental  and 
the  ENDyne  results  for  the  + H system.  Eneigy:  700  eV.  Basis  set:  pVDZ,  ENDyne 
transfer,  ’+’  ENDyne  elastic,  ’x’  transfer  and  ’o’  elastic  from  Ref.  [68]. 


Scattering  Angle,  0 (deg) 

Figure  4-14:  Reduced  differential  cross  sections  versus  the  scattering  angle.  Experimental  and 
the  ENDyne  results  for  the  H*  + H system.  Energy:  1000  eV.  Basis  set:  pVDZ.  ’*’  ENDyne 
transfer,  ’+’  ENDyne  elastic,  ’x’  transfer  and  ’o’  elastic  from  Ref.  [68]. 
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The  same  comments  and  conclusions  about  the  transfer  probability  hold  for  the 
reduced  differential  cross  sections.  We  should  only  add  that  the  agreement  between 
the  END-SD-FGWP  and  experimental  results  for  the  reduced  cross  sections  is  worse 
than  for  the  transfer  probability.  The  reduced  differential  cross  sections  are  explicitly 
dependent  upon  the  scattering  angle.  The  experimental  data  were  normalized  and  fitted 
to  theoretical  results,  which  are  based  upon  prescribed  trajectories.  We  can  then  argue 
that  the  experimental  results  for  the  reduced  cross  sections  have  been  downgraded  in 
quality  by  fitting  to  unsuitable  theoretical  results. 

The  H~^  + He  Collision 

The  behavior  of  proton-helium  scattering  is  different  from  that  of  the  proton-hydrogen 
collisions.  The  -f  He  system  exhibits  electron-electron  interaction,  the  charge  transfer 
is  non-resonant,  and  the  interaction  potential  has  repulsive  and  attractive  parts.  These 
facts  make  the  scattering  of  a proton  by  a helium  atom  an  important  system  to  study 
with  any  new  time-dependent  method  that  has  demonstrated  the  ability  to  describe  well 
the  -t-  H collision. 

Since  the  early  1930’s  the  proton-helium  scattering  has  been  studied  by  theoretical 
methods  [71].  The  accumulated  theoretical  data  falls  mostly  in  the  energy  range  above  10 
keV  and  are  related  mainly  to  total  (integral)  cross  sections.  Experimental  and  theoretical 
studies  of  differential  cross  sections  for  + He  are  more  scarce.  However,  it  is  possible 
to  find  several  experimental  studies  for  energies  down  to  5 keV  and  scattering  angles 
greater  than  1°  [72].  These  experiments  are  probing  the  repulsive  part  of  the  interaction 
potential  and  no  structures  are  then  observed.  Recent  studies  of  angular  differential 
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scattering  at  keV  energies  in  ion-atom  collisions  have  focused  on  very  small  angles, 
in  general,  below  1°  [73].  These  studies  are  motivated  by  the  highly  forward-peaked 
character  of  the  cross  sections  and  by  the  location  of  the  classical  rainbow  angle  within 
the  0-1°  range  (at  keV  energies).  The  fact  that  structures  occur  for  very  small  angle 
has  been  the  main  reason  for  the  delay  in  investigating  this  region,  since  high  angular 
resolution  is  needed.  From  the  theoretical  point  of  view,  the  oscillatory  patterns  seen  in 
elastic  and  inelastic  scattering  have  been  studied  extensively  using  semiclassical  scattering 
theories  [63,  74,  75]. 

Trajectory  Analysis 


Some  of  the  dynamical  trajectories  obtained  from  an  ENDyne  calculation  of  a proton 
being  scattered  by  a helium  atom,  with  collision  energy  of  50.0  eV  is  presented  in  Figure 
4-15.  In  general,  these  trajectories  should  be  expected  in  any  ion-atom  collision  non- 
resonant scattering. 
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Figure  4-15:  Dynamical  trajectories  of  being  scattered  by  He.  ENDyne  results  at  energy  of 
50.0  eV  with  a pVDZ  basis  set. 

The  appearance  of  classical  rainbow  scattering  implies  that  the  classical  deflection 
function  has  at  least  one  extremum.  This  means  that  there  are  at  least  two  impact 
parameters  that  will  yield  the  same  scattering  angle.  In  the  case  of  a proton  colliding 
with  a helium  atom  we  show  in  Figure  4-16  an  example  of  three  dynamical  trajectories 
that  yield  the  same  scattering  angle. 
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Figure  4-16:  Dynamical  trajectories  for  scattering  angle  of  2.65°.  ENDyne  results  for  + He 
system  at  50.0  eV  with  a pVDZ  basis  set.  Impact  parameters:  solid  line  ( — ) =1.12  a.u.,  dashed 
line  ( — ) = 1.60  a.u.,  dotted  line  (•  • •)  = 2.10  a.u. 

We  should  point  out  that  apparently  the  target  (He)  moves  only  in  the  perpendicular 
direction  (y-axis)  to  the  scattering  axis  (x-axis).  However,  this  is  due  to  the  scale  of  the 
graphic.  In  fact,  the  atomic  target  is  first  attracted  to  the  incoming  ionic  projectile  and 
then  repelled,  as  can  be  seen  in  Figure  4-17  for  three  selected  impact  parameters. 
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Figure  4-17:  Motion  of  the  atomic  target  (He)  during  a collision  with  at  50.0  eV  as  computed 
by  ENDyne  with  a pVDZ  basis  set.  The  initial  coordinate  of  the  He  atom  is  (15.0,  0.0).  Impact 
parameters:  solid  line  ( — ) = 1.12  a.u.,  dashed  line  ( — ) = 1.60  a.u.,  dotted  line  (•••)  = 2.10  a.u. 

The  Deflection  Function  and  the  Interaction  Potential 

In  general,  the  deflection  function  of  an  ion-atom  system  has  at  least  one  extreme 
point.  The  exception  is  the  resonant  collision  of  a proton  with  a hydrogen  atom.  This 
behavior  of  the  deflection  function  is  due  to  the  interaction  potential  which  for  the  non- 
resonant case  has  both  repulsive  and  attractive  parts.  We  present  a detailed  analysis  of 
the  deflection  function  for  the  + He  system  at  50.0  keV.  At  this  energy,  the  elastic 
scattering  is  the  only  important  process.  Consequently,  the  classical  elastic  scattering 
theory  can  be  used  to  invert  the  deflection  function  to  yield  the  intermolecular  potential. 
This  potential  is  compared  to  the  dynamic  potential  and  to  the  Bom-Oppenheimer 
potential  energy  surface.  However,  before  this  comparison  is  made,  we  show  the 
deflection  function  in  Figure  4-18. 
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Figure  4-18:  Deflection  function  for  + He.  ENDyne  results  at  50.0  eV  with  a pVDZ  basis  set. 

From  the  expression  for  the  classical  deflection  function,  equation  4,3,  it  is  possible 
to  devise  a procedure  that  takes  monotonic  parts  of  the  deflection  function  and  through  an 
inversion  yields  the  potential  function.  This  inversion  procedure  consists  of  the  following 
steps  [76,  63]: 


1.  define  some  distance  of  closest  approach  s. 


2.  evaluate  numerically  the  integral 


I{s)  = 7T 


1 


I 


^0 


0([s^  + 

\/ x'^ 


1 

/ 


y/y^S^  -I- 1/2 


(4.13) 


0 


3.  compute  r from 


r = 


(4.14) 
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4.  determine  a point  on  the  potential  energy  curve  (r,  V'(r))  from  the  above  expression 
for  r and  from  the  following  expression  for  y(r) 


V(r)  = E{1  - e"2^W). 


(4.15) 


5.  determine  values  of  V (r)  at  increasing  separations  r by  gradually  increasing  the  value 
of  the  lower  limit  s from  which  0(6)  is  integrated. 

In  order  for  this  inversion  procedure  to  be  valid,  the  following  conditions  must  be  met. 
(i)  Only  the  monotonic  branches  of  the  deflection  function  must  used.  In  the  case  of  a 
proton  colliding  with  helium,  the  deflection  function  must  be  divided  into  three  parts  to  be 
inverted,  namely,  {0, 6o),  {6q, K),  and  {6r, oo},  and  (ii)  E > K(ro),  that  is,  the  potential 
can  only  be  evaluated  up  to  the  classical  turning  point  ro  at  the  energy  E,  in  addition 
(iii)  only  data  that  do  not  involve  orbiting  collisions  may  be  inverted  uniquely,  that  is, 

,,/  X rdV(r) 


When  this  procedure  is  applied  to  the  deflection  function.  Figure  4-18,  generated  by 
ENDyne  for  + He  at  50.0  eV,  the  potential  curve  in  Figure  4-19  is  obtained. 
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Figure  4-19:  Potential  energy  curve  obtained  from  the  inversion  of  the  deflection  function 
generated  by  ENDyne  for  + He  at  50.0  eV  with  a pVDZ  basis  set. 

During  an  evolution  calculation  the  ENDyne  program  creates  a file,  called  restart  file, 
that  contains  the  history  of  the  system  evolving  in  time.  We  then  use  another  program, 
called  EVOLVE,  to  extract  relevant  information  about  the  evolution.  At  each  instant  in 
time  during  the  evolution  the  position  of  each  particle  and  the  potential  energy  of  the 
system  are  recorded.  As  a result,  we  can  compute  the  dynamical  potential,  that  is,  the 
potential  energy  at  each  instant  in  time  with  its  corresponding  interatomic  distance.  This 
dynamical  potential  can  be  obtained  from  data  before  and  after  the  collision,  and  also 
from  several  impact  parameters.  The  data  from  both  before  and  after  the  collision  of 
proton  -i-  helium  for  impact  parameters  ranging  from  0.01  to  1.0  a.u.,  yielded  the  same 
(up  to  six  decimal  places)  potential  curve.  One  of  this  curves  is  presented  in  Figure  4-20. 
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Figure  4-20:  Potential  energy  curve  obtained  from  the  direct  time  evolution.  ENDyne  results  for 
-I-  He  at  50.0  eV  with  a pVDZ  basis  set. 

We  also  used  an  ab  initio  CIS  (configuration  interaction  with  single  excitations)  to 
obtain  the  potential  energy  curves  for  the  ground  state  and  several  excited  states.  The 
results  for  + He  are  presented  in  Figure  4-21. 
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Figure  4-21:  Potential  energy  curves  for  + He  obtained  from  a CIS  calculation  with  a pVDZ 
basis  set. 

Attempts  to  plot  the  ground  state  of  a CIS  calculation,  together  with  the  inverted 
deflection  function  and  the  dynamical  potential  show  that  they  all  lie  on  top  of  each  other. 
This  is  a very  reassuring  result,  since  consistency  between  three  different  approach  has 
been  obtained.  The  reason  for  this  agreement  is  due  to  the  large  separation  between  the 
potential  curves  as  can  be  seen  in  Figure  4-21.  This  means  that  the  dynamics  will  remain 
on  only  one  surface  and  the  non-Bom-Oppenheimer  effects  are  negligible.  In  addition, 
unlike  the  H'*'  -h  H system,  the  energy  to  localize  the  electron  around  the  helium  atom 
in  the  incoming  asymptotic  channel  is  very  small  and  the  initial  state  basically  is  the 
ground  state  of  the  (H-t-He)"^  system. 
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Classical  Differential  Cross  Section 

Recent  improvements  in  data  acquisition,  analysis  and  detection  techniques,  error 
handling,  and  beam  collimation  have  enabled  experimentalist  to  explore  regions  of 
very  small  scattering  angle  [77,  78,  73].  The  angular  resolution  achieved  by  these 
improvements  is  now  good  enough  to  show  the  oscillatory  behavior  of  the  ion-atom 
differential  cross  sections.  Oscillations  in  the  differential  cross  sections  can  be  interpreted 
as  due  to  interferences  between  different  trajectories  leading  to  the  same  scattering 
angle.  Since  our  treatment  of  the  nuclear  motion  is  classical,  we  should  not  expect 
the  appearance  of  these  oscillations  in  the  classical  differential  cross  sections.  The  only 
structure  appearing  in  the  classical  treatment  is  due  to  the  oscillations  in  the  probabilities 
for  each  process. 

In  this  subsection,  we  present  the  theoretical  results  obtained  by  the  END-SD-FGWG 
method.  The  semiclassical  corrections  to  our  classical  dynamics  and  comparisons  with 
experimental  cross  sections  are  presented  in  the  next  subsection. 

We  present  the  classical  differential  cross  section  for  H'*’  + He  at  500  eV  in  some 
detail,  and  just  a summary  of  the  results  for  1500  and  5000  eV.  Since  the  differential 
cross  section  depends  upon  the  deflection  function,  we  start  by  showing  the  three  regions 
of  the  deflection  function  in  Figure  4-22.  The  regions  are, 

I)  {0, 6o } where  6q  corresponds  to  the  value  of  the  impact  parameter  where  the  deflection 
function  vanishes,  the  so-called  glory  angle,  0(6q), 

II)  ( 6o , 6r } where  W defines  the  rainbow  angle,  that  is,  the  extremum  of  the  deflection 


function 
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and, 

in)  {br,oo}  the  attractive  region. 


(4.17) 


Figure  4-22:  Deflection  function  for  + He.  ENDyne  results  at  500  eV  with  a pVDZ  basis 
set.  Regions:  I)  diamonds  ’o’,  II)  triangles  ’A’,  and.  III)  open  circles  ’o’. 

As  stated  before,  from  the  shape  of  the  deflection  function  in  Figure  4-22,  we  can 
see  that  more  than  one  impact  parameter  leads  to  the  same  scattering  angle  0.  Thus,  the 
expression  for  the  classical  differential  cross  section  becomes 


t 


bi  P{bi) 


sin  0 


dQ  I 
dh  l6=6i 


(4.18) 
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i.e.,  a sum  of  contributions  from  these  different  trajectories.  The  contributions  of  each  one 
of  the  three  regions  for  the  elastic  and  for  the  charge  transfer  differential  cross  sections 
are  shown  in  detail  in  Figures  4-23-4-24. 


Figure  4-23:  Elastic  differential  cross  section  for  H""  + He.  ENDyne  results  at  500  eV  with  a 
pVDZ  basis  set.  Regions:  I)  diamonds  ’o’,  II)  triangles  ’A’,  and,  III)  open  circles  ’o’. 
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Figure  4-24:  Charge  transfer  differential  cross  section  for  + He.  ENDyne  results  at  500  eV 
with  a pVDZ  basis  set.  Regions:  1)  diamonds  ’o’,  II)  triangles  ’A’,  and,  III)  open  circles  ’o’. 

The  elastic  and  charge  transfer  cross  sections  for  + He  at  500  eV,  1500  eV,  and 
5000  eV  are  shown  in  Figures  4-25-4-27. 
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Figure  4-25:  Elastic  and  charge  transfer  differential  cross  sections  for  + He.  ENDyne  results 
at  500  eV  with  a pVDZ  basis  set. 


Figure  4-26:  Elastic  and  charge  transfer  differential  cross  sections  for  H*  + He.  ENDyne  results 
at  1500  eV  with  a pVDZ  basis  set. 
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Figure  4-27:  Elastic  and  charge  transfer  differential  cross  sections  for  + He.  ENDyne  results 
at  5000  eV  with  a pVDZ  basis  set. 

As  we  can  see  from  Figures  4-25—4-27  the  classical  differential  cross  section  becomes 
singular  at  the  rainbow  angle,  and  tends  to  infinity  when  the  scattering  angle  goes  to  zero. 
We  also  should  notice  the  difference  in  orders  of  magnitude  between  the  elastic  and  the 
charge  transfer  processes. 

Semiclassical  Elastic  Differential  Cross  Section 


As  mentioned  before,  the  resolution  of  the  experimental  differential  cross  sections 
is  high  enough  to  allow  us  to  observe  quantum  effects.  Consequently,  the  results 
presented  for  the  classical  differential  cross  sections  cannot  directly  be  compared  with 
the  experimental  data,  except  for  the  value  of  the  rainbow  angle.  In  order  to  remedy 
this  problem  we  use  semiclassical  corrections  [63,  75,  74,  79,  80]  to  our  classical  elastic 


differential  cross  sections. 
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We  start  by  reviewing  one  of  the  old  stablished  techniques  to  solve  the  elastic 
scattering  problem  using  time-independent  methods  . The  elastic  scattering  by  a spherical 
potential  V{r)  can  be  described  by  the  wave  function  ^(f)  which  is  an  eigenfunction 


of  the  Schrodinger  equation 


_ ^y(r)  -f  ^{f)  = 0 


(4.19) 


with 


p = hk,  k^  = 


(4.20) 


where  p is  the  momentum.  The  elastic  differential  cross  section  can  be  expressed  as 


follows 

^ = \Am 

where  A{0)  is  called  the  scattering  amplitude. 


(4.21) 


In  order  to  obtain  the  scattering  amplitude  A{0)  we  expand  the  solution  of  the 


SchrOdinger  equation  for  the  scattering  problem  in  terms  of  Legendre  polynomials 


. OO 

^(^  = - Citpi{r)Pi{cos  6) 


(4.22) 


/=o 


where  C/  are  the  amplitudes  of  the  partial  waves  ipi{r),  Pi{cos0)  are  the  Legendre 
polynomials,  and  / is  the  angular  momentum  quantum  number.  This  is  the  so  called 
partial  wave  expansion,  and  it  yields  the  following  set  of  radial  equations  for  V’/(r) 


- ?^V(r)  - 
h H 


V’/(r)  = 0 . 


(4.23) 


Using  the  normalization  condition  for  ^'(f)  and  the  orthogonality  of  the  Legendre 
polynomials,  we  find  that 

..  OO 

^ /=0 


(4.24) 
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where  the  phase  shift  is  defined  as 


6i  = lim  {Phase  tpi{r)  — Phase  ‘4’i{r,V{r)  = 0)} 


(4.25) 


Consequently,  the  crux  of  quantum  mechanical  scattering  the  calculation  of  the  phase 
shift.  Although  the  expression  for  the  scattering  amplitude,  equation  (4.24),  is  exact,  the 
series  converge  extremely  slowly.  Several  approximate  techniques  have  been  developed 
that  give  analytical  expressions  for  the  scattering  amplitude. 

For  instance,  employing  the  WKB  semiclassical  approximation,  the  following  ex- 
pression for  the  phase  shift  can  be  obtained  [75], 


OC 


(4.26) 


where 


(4.27) 


and 


(4.28) 


Differentiating  the  WKB  phase  shift  6i  with  respect  to  / we  obtain 


OO 


(4.29) 


To 

which  when  compared  with  the  classical  deflection  function  (a  = tq) 


OO 


OO 


2J  / drr-^FI~^^^{r) 


(4.30) 
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suggests  that 


(4.31) 


This  shows  that  the  semiclassical  phase  shift  is  connected  to  the  classical  deflection 
function  of  the  same  angular  momentum.  Now  that  we  have  an  analytical  expression 
for  the  semiclassical  phase  shift  we  can  use  approximations  of  the  Legendre  polynomials 
Pi{cos0)  to  obtain  expressions  for  the  scattering  amplitude.  We  analyze  three  cases, 
namely,  a)  angles  not  too  small  and  not  close  to  the  rainbow  angle  Or,  b)  very  small 
angles,  and  c)  angles  very  near  to  Or. 

(a)  {kb)~^  < 0 < TT {4  and  0 / Or- 

This  range  covers  scattering  at  most  angles,  excluding  very  small  angles  and  near  the 
rainbow  angle  Or.  In  this  range,  the  Legendre  polynomials  P/(cos  0)  are  approximated  by 


2 


cos [(/ -t- 1/2)0  — 7r/4]  + 0(/  ^) 


(4.32) 


and  the  WKB  scattering  amplitude  can  be  simplified  to 


(4.33) 


0 


with 


(4.34) 


Because  of  the  small  value  of  h these  integrals  may  be  evaluated  by  the  method  of 
stationary  phase.  After  some  algebraic  manipulations,  the  semiclassical  (sc)  elastic 
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differential  cross  section  is  given  by 


-SC  

= 


da 


2 cos(a,  - oy)  + E(‘'5')i 


(4.35) 


• • 
X<] 


where 


(4)i  = - 


bi 


sin^ 


d© 

db 


(4.36) 


b=bi 


is  the  classical  cross  sections  originating  from  each  impact  parameter  that  gives  rise  to 
the  same  scattering  angle.  Assuming  that  -tt  < 0 < x,  the  phase  difference  can  be 
expressed  as 


Aq^  = a,-  — Qj 


bi 


= k J db  0(6)  - ke{r]ibi  - rjjbj)  + ^ [(77, • - 27/^;)  - {T)j  - 2t)s.)] 
where  rjs  is  the  algebraic  sign  of  the  slope  of  the  deflection  function 


(4.37) 


(4.38) 


and 

, + 1 if  0(6,)  > 0 

^ . (4.39) 

^ - 1 if  0(6,)  < 0 

We  find,  therefore,  that  unless  the  deflection  function  0 is  a monotonic  function  of  the 
impact  parameter  6 and  is  confined  to  values  between  — x < 0 < x (so  that  for  a 
given  6 there  is  but  one  value  of  6),  we  do  not  get  the  classical  result  for  the  cross 
section,  no  matter  how  short  the  wavelength.  Instead  there  are  additional  interference 
terms.  However,  for  truly  macroscopic  scattering  the  beam  of  projectiles  is  not  exactly 
monoenergertic.  Thus,  since  the  phase  angles  a,-  contain  the  microscopic  variable  k,  or 
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the  wavelength,  a spread  in  the  energy,  gives  averaged  phase  angles  yielding  the  purely 
classical  result. 

(b)  Very  small  angles:  9 < {kb)~^  ~ 0 and  9 ^ 9^. 

This  range  of  scattering  angle  is  called  forward  glory  scattering.  One  of  the  classical 
trajectories  which  contributes  to  the  forward  glory  corresponds  to  a large  impact  parameter 
(6i  in  Figure  4-18).  The  other  is  the  glory  trajectory  at  small  b,  corresponding  the 
coalescence  of  62  and  63  trajectories  in  Figure  4-18.  The  asymptotic  expression  for  the 
Legendre  polynomials  Pi{cos9)  at  very  small  angles  is 


0 

As  a result,  the  scattering  amplitude  for  very  small  angles  can  be  expressed  as 


P/(cos  9)  ~ (cos  9yj()  [(/  -|- 1 /2)  sin  9\  ~ (cos  9)^Jq  [(/  -|- 1 /2)0]  (4.40) 


where  the  Bessel’s  integral  is 


(4.41) 


(4.42) 


We  expand  about  the  glory  value,  (0(6o)  = 0, 


e{i)  = e{io)  + -^  , , (/-/o)  + ... 

ai  /=/o 


(4.43) 


and  obtain  the  scattering  amplitude 


db  6=60 


(4.44) 
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TTius,  the  elastic  differential  cross  section  for  angles  close  to  the  forward  glory  angle 
is  given  by 


de 


dh 


b=bo 


[Jo(^^osin  0)]^. 


(4.45) 


(c)  Near  the  rainbow  angle:  B k Or- 


The  deflection  function  has  an  extremum  at  0r,  so  we  expand  0 about  it 


0 = 0r  + 


1 SQ 


db^ 


b=br 


so  that. 


c c 1 ^ / T T ^ 1 d^0 

h ^ h, + 56,(7  - 7,)  + 


b=br 


(J-Jr)' 


The  scattering  amplitude  close  to  the  rainbow  angle  is  given  by 


A{B)  ~ , 


1 


2irkb 
sin  $ 


nl/2 


,«[a+»/fc6r(Sr-^)] 


Ai(.) 


where 


7T 


Q = 26b  - kbrQr  + TJT 

4 


and  the  Airy  function  is  defined  as 


00 


Ai(2)  = — J ducos{zu  + 


0 


with  the  argument  given  by 


2 = 2^l^k^l^1)rn{0T  - B) 


<Pe 


-1/3 


b=br 


Vt  = sgn 


<i62 


b=br 


(4.46) 


(4.47) 


(4.48) 


(4.49) 


(4.50) 


(4.51) 
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The  value  of  the  elastic  differential  cross  section  near  the  rainbow  angle  can  then  be 
written  in  the  form 


25/Vibi/36r 

sin  6 


(Pe 

d9 


-2/3 


|Ai(2))^ 


(4.52) 


The  expression  in  equation  (4.52)  shows  the  complex  structure  of  the  elastic  differential 
cross  section  in  the  region  of  Br  when  compared  to  the  classical  differential  cross  section. 
That  is,  the  classical  singularity  at  ^ has  been  smoothed  out  and  replaced  by  an 
exponential  damping  2tX  B > Or  and  by  an  oscillatory  feature  at  ^ This  oscillating 
behavior  arises  from  the  Airy  function  for  ^ < 0 and  the  maxima  in  the  differential  cross 
section,  called  rainbow  maxima,  correspond  to  the  extrema  in  Ai(2),  that  is,  2 = -1.0188, 
-3.2482,  -4.8201,  -6.1633,  etc.  Since,  in  general,  0(6  « 6r)  < 0 and  d‘^Q/db'^\i,=i^  > 0 


we  have  that 


2 = -{Br  - B) 


1 d'^Q 


2P  d9 


-1/3 


h=hr 


= {B-  Br)q-^l^ 


(4.53) 


1 <PQ 


9 = 


2it2  dip 


b=br 


Thus,  the  main  peak  occurs  below  Br,  where  2 = -1.0188,  i.e. 


Omax  =Br-  I.OISS?^/^  (4  54) 

Subsidiary  supernumerary  rainbows  occur  at  lower  angles  (2  = -3.2482,  -4.8201, 
-6.1633,  etc). 

Experimental  and  Theoretical  Elastic  Differential  Cross  Sections 


The  value  of  the  rainbow  angle,  Br,  can  be  determined  experimentally  as  the  position 
of  the  inflection  in  Od{B)  beyond  the  rainbow  maximum.  For  the  scattering  of  a proton 
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with  a helium  atom  the  experimental  determination  of  the  rainbow  angle  is  compared 
with  the  theoretical  result  in  Table  4-5.  The  agreement  between  the  ENDyne  results  and 
the  experimental  rainbow  angles  is  excellent. 


Table  4-5:  Experimental  and  theoretical  rainbow  angles  for  + He  system. 


Energy 

(eV) 

Rainbow  angle  (degrees) 

Impact  parameter 

(a.u.) 

ENDyne 

Experimental® 

ENDyne 

50.0 

— 

2.963 

1.826 

500 

0.32 

0.3015 

1.778 

1500 

0.11 

0.1013 

1.772 

5000 

0.03 

0.0302 

1.772 

a:  Ref.  [73]. 


We  are  working  in  the  angular  range  of  0.01° — 1.0°  and  0.04°  < {kb)~^  < 0.29° 
for  H'*’  + He  at  500  eV.  Consequently,  we  are  in  the  regime  of  very  small  angles, 
^ ~ 0 and  should  use  equation  (4.45)  for  the  differential  cross  section  far  from 

the  rainbow  angle.  The  semiclassical  corrections  to  the  classical  elastic  differential  cross 
section  obtained  by  the  END-SD-FGWP  method  is  shown  in  Figure  4-28. 
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Figure  4-28:  Elastic  differential  cross  section  for  + He.  Energy:  5000  eV.  Basis  set:  pVDZ. 
Solid  line  corresponds  to  the  semiclassical  corrections  to  the  ENDyne  results.  Open  circles  are 
experimental  data  [73]. 

We  have  presented  the  results  for  500  eV  only  since  at  this  energy  the  elastic  process 
is  several  orders  of  magnitude  larger  that  the  charge  transfer.  Thus,  the  semiclassical 
corrections  to  the  classical  elastic  differential  cross  section  should  be  good,  as  can  be 
seen  from  the  comparison  between  the  theoretical  and  the  experimental  results.  It  remains 
to  develop  an  appropriate  formalism  within  END  that  yield  these  quantum  interferences 
for  inelastic  processes. 

The  H2  Collision 

We  have  seen  so  far,  that  the  END  formalism  based  upon  a single  determinant  and 
classical  nuclei  gives  very  good  results  for  resonant  and  non-resonant  ion-atom  scattering. 
We  now  examine  the  performance  of  this  method  for  ion-molecule  collisions,  in  particular 

+ H2.  In  these  type  of  collisions  not  only  do  the  electrons  play  an  important  role 
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(transfer  and  excitation)  but  so  do  the  nuclei,  since  vibrational  and  rotational  excitations 
take  place.  There  is  a surprising  lack  (virtually  none)  of  rigorous  theoretical  investigations 
on  dynamics  which  occur  in  ion-molecule  collisions  in  low-  to  intermediate-energy 
regime  (50  eW  < E < 50  keV).  There  are  several  reason  for  this  lack  of  theoretical 
treatment,  but  mainly  that,  it  is  quite  a complex  problem  to  obtain  reasonably  accurate 
adiabatic  potentials  and  eigenstates  as  functions  of  intemuclear  distances  and  molecular 
orientations.  In  addition,  for  an  ion-molecule  system,  the  number  of  internal  degrees 
of  freedom  that  need  a proper  dynamical  treatment  increases  dramatically.  The  END 
formalism  in  this  context  gains  a lot  of  importance  since  it  does  not  need  adiabatic 
potentials  and  treats  properly  all  degrees  of  freedom.  In  the  present  formulation,  namely, 
END-SD-FGWP,  we  have  not  incorporated  the  ETF’s  (electron  translation  factors),  so 
it  is  limited  to  low-energies  (E  < I keV)  applications.  The  END  formalism  is  also 
capable  of  treating  dynamics  at  very  low-energies  (E  < 10  eV),  which  will  be  presented 
in  Chapter  5. 

In  this  subsection  we  present  the  collision  of  a proton  with  a hydrogen  molecule  in 
its  electronic  and  vibrational  ground  state  at  500  eV.  We  examine  some  of  the  extreme 
molecular  orientations,  namely,  (0°,  0°)  and  (90°,  0°).  This  should  be  representative  of 
the  effects  of  orientation  of  the  target  molecule  on  the  charge-transfer  mechanism,  cross 
sections  (charge  transfer),  and  energy  transfer  mechanisms. 

Trajectory  Analysis 

The  nuclear  motion  of  an  ion-molecule  collision  is  quite  a bit  more  complicated  (and 
more  interesting)  than  the  ion-atom  scattering  showed  in  the  last  two  subsections.  In  an 


89 


ion-molecule  collision  one  would  expect  that  vibrational  and  rotational  excitations  occur. 
These  excitations  are  shown  graphically  via  the  evolution  of  the  nuclear  coordinates 
in  time.  Figure  4-29  shows  the  head-on  collision  (0°,  0°  orientation)  for  three  impact 
parameters  that  lead  to  the  same  scattering  angle  (0.29°).  We  should  clarify  the  notation 
for  the  molecular  orientation  (a,  /3).  The  first  angle,  a,  is  the  angle  between  the  scattering 
direction  and  the  molecule  axis,  and  is  the  dihedral  angle  between  the  scattering  plane 
and  the  molecular  plane.  Figure  4-29  represents  the  nuclear  motions  in  the  laboratory 
frame  where  the  proton  trajectories  are  the  lines  (almost  straight  lines)  and  the  motion 
of  the  hydrogen  atoms  in  the  H2  molecule  are  the  oscillatory  lines.  A more  detail 
representation  of  these  motions,  namely,  the  target  atom  coordinates,  are  plotted  in 
Figure  4-30. 


Figure  4-29:  Dynamical  trajectories  for  + H2  collision  from  ENDyne  calculation  for  (0°,  0°) 
orientation  at  500  eV  with  a pVDZ  basis  set. 
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Figure  4-30:  Dynamical  trajectories  of  the  taiget  H2.  ENDyne  results  for  (0°,  0°)  orientation  at 
500  eV  with  a pVDZ  basis  set. 

The  nuclear  motion  in  the  laboratory  frame  for  the  (90°,  0°)  orientation  is  plotted  in 
Figure  4-31.  These  trajectories  lead  to  the  same  scattering  angle  (0.29°)  as  in  the  case 
of  (0°,  0°)  orientation.  As  we  can  see,  the  trajectories  for  both  molecular  orientations 
are  quite  distinct,  but  a closer  look  shows  the  same  qualitative  behavior,  i.e.,  at  small 
impact  parameters  the  target  and  projectile  initially  attract  each  other  then  the  repulsive 
interaction  overcomes  this  attraction  and  they  move  apart  from  each  other.  At  large 
impact  parameters  they  attract  each  other  and  move  close  together,  but  the  target  cannot 
keep  up  with  the  projectile  velocity.  In  both  cases  the  H2  molecule  starts  vibrating  and 
rotating  after  the  collision.  We  should  notice  that  even  for  the  longest  evolution  times, 
namely,  1200  a.u.  (=  30  fs)  for  (0.0°,  0.0°)  and  800  a.u.  (=  20  fs)  for  (90.0°,  0.0°), 
considered  here,  the  rotational  motion  is  almost  absent. 
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X-Coordinate  (a.u.) 


Figure  4-31:  Dynamical  trajectories  for  H"^  + H2  collision  from  ENDyne  calculation  for  (90°, 
0°)  orientation  at  500  eV  with  a pVDZ  basis  set. 


X'Coordinate  (a.u.) 


Figure  4-32:  Dynamical  trajectories  of  the  target  H2.  ENDyne  results  for  (90°,  0°)  orientation 
at  5(X)  eV  with  a pVDZ  basis  set. 
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Vibrational  Analysis 


The  trajectory  analysis  seems  to  show  that  the  H2  molecular  target  vibrates  in  quite 
different  vibrational  states  (frequency  and  amplitude).  We  emphasize  again  that  the 
trajectories  shown  in  the  last  subsection  lead  to  the  same  scattering  angle  (0.29°).  The 
angular  dependence  of  the  vibrational  populations  p„  of  H2  molecules  after  a collision 
with  H''^  have  been  extensively  studied  by  theoretical  and  experimental  methods  for  low- 
(E  < 30  eV)  and  high-energies  (E  > 10  keV).  However,  due  to  recent  experimental 
developments  vibrationally  resolved  measurements  [81],  differential  in  scattering  angle 
6,  have  been  made  for  the  intermediate-energy  range  (100  eV  < < 10  keV)  and  started 

challenging  the  theoreticians.  Since  we  have  a classical  description  of  the  nuclear  motion, 
some  approximate  vibrational  wavefunction  would  have  to  be  determined  to  yield  p„  as  a 
function  of  6.  We  have  not  fully  analyzed  this  problem  yet,  so  we  present  in  this  section 
the  vibrational  analysis  of  the  classical  motion  of  the  nuclei. 

We  find  that  the  best  way  to  obtain  vibrational  frequencies  of  the  H2  target  is  to  use  the 
Prony  method  [82]  to  analyze  the  interatomic  distance  of  the  H2  molecule  as  function  of 
time.  Figure  4-33  shows  the  interatomic  distances  for  all  six  impact  parameters  that  lead 
to  the  same  scattering  angle  (0.29°)  for  both  (0°,  0°)  and  (90°,  0°)  molecular  orientations. 
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Figure  4-33:  Interatomic  distances  of  the  target  H2  as  function  of  time.  ENDyne  results  for  (0°, 
0°)  and  (90°,  0°)  orientations  at  5(X)  eV  with  a pVDZ  basis  set.  (1  a.u.  = 0.024195  fs). 

The  Prony  method  is  an  alternative  to  the  DFT  (discrete  Fourier  transform)  for  spectral 
analysis  and  is  widely  used  in  signal  processing.  It  is  a technique  for  modeling  sampled 
data  as  linear  combination  of  (complex)  exponentials  with  damping  coefficients  e and 
angular  frequencies  uj 

V 

i(n)  = (4.55) 

i=\ 


where  the  number  of  terms  p is  called  the  order  of  the  method  and  the  equally  spaced 
steps  of  the  data  sample  is  n{6t).  Details  of  the  algorithm  and  its  implementation  can  be 
found  in  reference  [82].  Once  the  parameters  of  the  exponential  model  are  determined, 
the  frequency  spectrum  is  obtained  directly,  using  the  infinite  domain  analytic  Fourier 
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transform 


N-l 

x{k)  = N{6t) 

n=0 


(4.56) 


where  N is  the  number  of  points  in  the  data  sample.  In  order  to  enhance  the  signal  it 
is  useful  to  use  the  relative  spectral  densities  (RSD)  as  a function  of  the  frequency  u to 
obtain  the  spectrum.  The  RSD,  expressed  in  decibel  (dB),  is  defined  as 

RSDHldbl  = 101og.„(J^).  (4.57) 


The  vibration  spectrum  of  H2  for  a specific  trajectory  (b  = 1.25  a.u.  and  (a,  = 

(0°,  0°))  is  shown  in  Figure  4-34.  The  other  trajectories  present  similar  spectra.  The 
largest  order  p of  the  method  before  linear  dependency  occurs  was  around  40  for  all 
trajectories  analyzed. 


Figure  4-34;  Vibrational  spectrum  of  the  target  H2.  ENDyne  results  for  (0°,  0°)  orientation  with 
impact  parameter  of  1.25  a.u.  Eneigy:  500  eV.  Basis  set:  pVDZ. 
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In  Table  4-6  we  present  the  vibrational  frequencies  for  all  trajectories  shown  in  Figure 
4-33  (impact  parameters  leading  to  the  same  scattering  angle  (0.29°)  for  both  (0°,  0°)  and 
(90°,  0°)  molecular  orientations).  The  differences  in  the  vibrational  states  are  basically 
due  to  differences  in  the  energy  transfer  and  the  electronic  state  of  the  H2  target  after 
the  collision. 


Table  4-6:  Vibrational  frequencies  (cm  *)  as  function  of  the  impact  parameter  (a.u.)  and 
molecular  orientation  (a°,  13°). 


b (a.u.) 

0 

0 

Vibrational  Frequency 
(cm'^) 

1.25 

0,0 

4441 

1.59 

0,0 

4467 

3.42 

0,0 

4612 

1.70 

90,0 

4415 

2.20 

90,0 

4308 

3.50 

90,0 

4580 

Transfer  Probability 


Since  the  vibrational  ftequencies  are  different  for  a given  scattering  angle,  it  suggests 
that  the  electronic  states  of  the  H2  target  are  also  different.  The  transfer  probability 
reflects  these  differences  and  is  plotted  in  Figure  4-35. 
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Impact  Parameter  (a.u.) 


Figure  4-35:  Probability  for  charge  transfer  in  + H2  collision.  ENDyne  results  for  (0°,  0°) 
and  (90°,  0°)  orientations.  Eneigy:  500  eV.  Basis  set:  pVDZ. 

It  is  interesting  to  note  that  even  in  these  two  extreme  molecular  orientations  the 
transfer  probability  has  similar  qualitative  behavior,  however  the  quantitative  differences 
are  large  enough  to  yield  a variety  of  vibrational  frequencies  for  the  same  scattering  angle 

Deflection  Function 

As  we  have  seen  before,  the  deflection  function  is  an  indicator  of  the  appearance 
of  the  intermolecular  potential.  The  adiabatic  potential  curves  as  a function  of  the 
- H2  (center  of  mass)  distance  are  highly  dependent  upon  the  molecular  orientation  for 
intermolecular  distances  smaller  that  2.5  a.u.  [83].  We  should  expect  then  some  signif- 
icant differences  between  the  deflection  function  for  (0°,  0°)  and  (90°,  0°)  orientations 
at  small  impact  parameters.  These  deflection  functions  are  shown  in  Figure  4-36.  They 
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do  confirm  the  findings  of  ab  initio  calculations  for  the  dependence  of  the  potential  on 
the  molecular  orientation. 


Figure  4-36:  Deflection  functions  for  + H2  collision.  ENDyne  results  for  (0°,  0°)  and  (90°, 
0°)  orientations.  Energy:  500  eV.  Basis  set:  pVDZ. 

Differential  Ooss  Sections 

Since  we  have  the  probabilities  and  the  deflection  function  we  can  compute  the 
classical  differential  cross  section  for  these  two  molecular  orientation.  The  calculation  of 
the  differential  cross  sections  is  also  motivated  by  the  experimental  measurements  [84] 
of  the  transfer  and  elastic  processes  in  + H2  collision. 

Figures  4-37—4-38  show  the  elastic  and  charge  transfer  differential  cross  sections 
calculated  by  ENDyne  for  the  two  molecular  orientations  considered  here  and  compared 
with  the  experimental  data  [84]. 
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Figure  4-37:  Elastic  differential  cross  section  for  + H2  collision.  ENDyne  results  for  (0°, 
0°)  and  (90°,  0°)  orientations.  Experimental  data  from  reference  [84].  Energy:  500  eV.  Basis 
set:  pVDZ. 


Figure  4-38:  Charge  transfer  differential  cross  section  for  + H2  collision.  ENDyne  results 
for  (0°,  0°)  and  (90°,  0°)  orientations.  Experimental  data  from  reference  [84].  Energy:  500 
eV.  Basis  set:  pVDZ. 
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The  calculated  classical  elastic  differential  cross  sections  for  both  molecular  orien- 
tations agree  quite  well  with  the  experimental  results.  From  the  inflection  in  the  exper- 
imental data  for  the  elastic  cross  section  we  estimate  that  the  rainbow  angle  should  lie 
between  4.0°  and  4.7°,  which  agrees  quite  well  with  the  ’averaged’  theoretical  rainbow 
angle  from  the  molecular  orientations  considered  here.  The  agreement  with  the  experi- 
mental data  for  the  charge  transfer  cross  section  is  not  as  good  as  for  the  elastic  one.  The 
reason  perhaps  is  that  we  are  comparing  classical  cross  sections  with  the  experimental 
results.  The  latter,  by  nature,  includes  the  quantum  effects  due  to  interference  between 
different  trajectories  leading  to  the  same  scattering  angle.  In  the  case  of  a molecular 
target,  these  effects  should  be  even  greater  than  for  atomic  target,  since  we  have  also 
interference  between  trajectories  from  different  molecular  orientations. 

Summary  and  Conclusions 


Both  quantum  and  classical  dynamical  quantities  calculated  by  the  END-SD-FGWP 
method  is  very  satisfactory.  Some  of  the  quantum  dynamical  quantities  like  transfer 
probability  and  interparticle  interactions  are  compared  with  experimental  results  and  it 
shows  an  excellent  agreement.  So  do  some  classical  quantities  like  rainbow  angles. 
However,  some  work  has  to  be  done  to  develop  appropriate  techniques  and  formalisms 
to  correct  classical  results  like  vibrational  populations  and  quantum  effects  in  differential 
cross  sections. 

The  systems  treated  here  + H,  He,  and  H2  present  a wide  variety  of  situations. 
For  instance,  we  have  resonant  (H'*'  + H),  near-resonant  (H"^  -1-  H2),  and  non-resonant  (H'*’ 
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+ He)  situations,  and  the  END-SD-FGWP  method  represented  them  quite  well.  Also, 
we  have  shown  that  the  molecular  orientation  plays  an  important  role  in  charge  transfer 
processes  as  do  the  quantum  effects  due  to  contributions  from  different  trajectories. 


CHAPTER  5 

INTRAMOLECULAR  CHARGE  TRANSFER  DYNAMICS 


The  dynamics  of  inter-  and  intramolecular  charge  transfer  has  been  a very  important 
subject  in  chemical-physics  sciences  [85-89].  The  main  reason  is  that  charge  transfer  is 
the  essence  of  some  most  fundamental  biological  and  biochemical  processes.  Also,  charge 
transfer  is  important  in  many  fields  such  as  atomic  physics,  plasma  physics,  astrophysics, 
semiconductor  physics,  organic  and  inorganic  chemistry,  and  many  others. 

In  this  chapter  we  study  the  following  charge  transfer  systems  Li — H-Li'*'  ^ Li"^- 
H — Li  and  Li-CN-Li"^  ^ Li'^-CN-Li.  These  molecular  systems  represent  a good  prototype 
for  theoretical  studies  of  charge  transfer  dynamics.  They  are  small  enough  to  be  treated  at 
some  high  level  ab  initio  time-dependent  and  time-independent  methods.  In  addition,  the 
Li-CN-Li  system  can  be  seen  as  exhibiting  charge  transfer  between  two  mixed  valence 
metal  atoms  mediated  by  a bridge. 

Electron  Transfer  Formalisms 

Due  to  its  conceptual  simplicity  the  formalism  for  electron  transfer  proposed  by 
Marcus  [90-92,  89]  in  the  late  50’s  has  drawn  a lot  of  attention  from  both  theoreticians 
and  experimentalists.  This  theory  has  been  extensively  reviewed,  revised,  and  extended. 
However,  two  basic  assumptions  are  still  being  employed,  namely,  a)  that  there  is  a 
reaction  coordinate  that  takes  the  reactants  to  the  products  via  nuclear  motion,  and  b) 
there  is  a coupling  (//12)  between  two  electronic  states  (donor  and  acceptor).  A typical 
electron  transfer  reaction  coordinate  is  shown  in  Figure  5-39. 
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Figure  5-39:  Diabatic  and  adiabatic  potential  curves  for  normal  electron  transfer. 

Non-dynamical  type  approaches  for  electron  transfer  assume  that  the  transfer  rate  k 
is  usually  given  by  an  expression  based  on  the  Fermi  golden  rule  [92,  89,  93] 

2ir 

»:  = y|/ri2p(FC).  (5.1) 

The  quantity  (FC)  is  the  Franck-Condon  factor  and  is  related  to  the  vibrational  spectrum 
of  the  donor-acceptor  system  and  its  surroundings,  if  any.  Analytical  expressions  for  (FC) 
derived  from  classical  and  semiclassical  theories  exist  and  are  discussed  in  references  [92, 
89].  The  matrix  element  H\2  represents  the  electronic  coupling  of  the  donor  and  acceptor 
(or  reactant  and  product).  This  is  the  quantity  that  has  been  challenging  theoreticians  and 
experimentalists.  There  are  several  approaches  to  compute  the  electron  transfer  matrix 
element.  Usually,  for  large  systems  such  as  proteins,  the  transfer  dynamics  is  treated  as 
a one-electron  problem  [87,  93].  However,  for  small  molecular  systems  many-electrons 
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ab  initio  and  semiempirical  theories  have  been  used  [89].  In  this  case,  the  most  common 
approach  to  calculate  H\2  is  to  compute  the  diabatic  states  which,  in  general,  are  non- 
orthogonal  and  use  them  as  zeroth  order  states  in  a perturbational  approach  [94,  89]  or 
use  them  in  a two-dimensional  non-orthogonal  Cl  (configuration  interaction)  problem 
[95,  96].  More  specifically,  the  broken  space  symmetry  wave  functions,  |1)  and  |2),  are 
used  in  the  following  expression  for  the  coupling  element  [97,  94,  96,  89] 


1 


1 Si2 


{\\H\2)-Sn 


(l|/f|l)-H(2|i/l2) 


(5.2) 


5i2  = (1|2). 


Since  the  nuclear  motion  is  responsible  for  electron  transfer  [97,  89],  some  attempts 
have  been  made  to  treat  the  motion  of  the  nuclei  explicitly.  Most  of  these  dynamical 
approaches  are  also  two-state  models  with  a coupling  H\2.  One  of  these  approaches  uses 
an  ab  initio  method  to  calculate  the  coupling  H12  at  each  nuclear  configuration  as  the 
system  is  evolving  in  time  [98,  99].  Some  simpler  approaches  use  a model  potential 
function  for  the  electronic  state  and  the  nuclear  motion  is  treated  explicitly  by  the  time- 
dependent  variational  principle  [100]  or  by  the  semiclassical  approach  of  Nikitin  [101]. 
All  these  approaches  for  electron  transfer  need  to  compute  the  coupling  matrix  element. 
In  order  to  obtain  H12  it  is  necessary  to  know  first  what  is  the  reaction  coordinate,  and 
that  might  be  very  difficult  to  determine  for  large  systems. 

We  use  the  ENDyne  program  based  on  the  END-SD-FGWP  formalism  to  study 
intramolecular  electron  transfer.  The  main  advantage  is  that  we  do  not  need  to  compute 
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the  coupling  matrix  element,  and  to  initialize  the  electron  transfer  we  just  distort  the 
molecule  and  let  the  system  evolve  in  time. 

Structure,  Energetics,  and  Electron  Density  of  LiCNLi 


Before  presenting  the  results  for  LiCNLi  we  should  specify  that  all  results  presented 
in  this  subsection  were  obtained  using  the  3-21G  basis  set  [102].  In  addition,  we 
perform  time-independent  ab  initio  self-consistent  field  (SCF)  calculations  to  determine 
potential  surfaces  and  electronic  densities  of  the  LiCNLi  and  LiHLi  molecules.  We 
do  not  describe  the  details  of  the  Hartree-Fock  method  and  the  methods  that  introduce 
electronic  correlation,  since  we  use  standard  approaches  widely  tested  and  reviewed  in 
the  literature  [103-105]. 

The  potential  energy  surface  of  the  LiCNLi  molecule  has  been  recently  studied  where 
the  basis  set  and  electronic  correlation  effects  were  studied  in  some  detail  [106].  We 
augmented  those  results  by  adding  data  from  a UHF/3-21G  calculation.  The  linear 
structure  of  LiCNLi  has  two  minima  at  all  levels  of  theory  explored.  Table  5-7  presents 
the  results  for  these  two  linear  structures  at  several  levels  of  theory.  As  we  can  see 
the  general  geometrical  structure  of  linear  LiCNLi  is  quite  independent  of  basis  set  and 
electronic  correlation  corrections. 
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Table  5-7:  Minima  for  the  linear  structure  of  LiCNLi  at  various  levels  of  theory.  Bond  distance 
in  pm. 


Theory  Level 

Structure  P 

Structure  IF 

1 

LiC 

CN 

NLi 

LiC 

CN 

NLi 

UHF/3-21G 

198.0 

115.0 

192.6 

214.1 

115.8 

178.8 

UHF/6-31G* 

198.3 

114.5 

198.8 

216.6 

115.3 

182.1 

UHF/6-31+G* 

196.5 

114.6 

197.9 

215.9 

115.4 

180.7 

UMP2/6-31G* 

196.0 

118.4 

200.1 

213.1 

118.8 

183.6 

CASS/3-21G 

197.0 

117.4 

194.3 

212.0 

117.7 

180.5 

a:  All  results  from  Ref.  [106]  except  UHF/3-21G. 


These  two  linear  structures  are  local  minima  in  the  potential  energy  surface.  The 
proposed  global  minimum  has  a C2v  structure  at  UHF/3-21G  level  with  the  lithium 
atoms  bonded  to  the  nitrogen.  Another  interesting  fact  about  the  potential  energy  surface 
of  this  system  is  that  the  transition  state  between  the  linear  structures  (I  and  II)  has  Cg 
symmetry.  The  structure  of  this  transition  state  and  global  minimum  at  the  UHF/3-21G 
level  is  presented  in  Figure  5-40. 


Transition  Slate  I-D 


Figure  5-40:  Transition  state  for  structures  I and  II  and  global  minimum  structures.  Bond 
distances  in  pm  and  angles  in  degrees. 


The  energetics  of  the  extrema  of  the  potential  energy  surface  of  LiCNLi  is  indicated 
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in  Table  5-8.  As  we  can  see,  the  qualitative  and  quantitative  features  of  the  potential 
energy  surface  of  LiCNLi  are  quite  basis  set  independent,  and  only  the  relative  stability 
of  the  linear  structures  depends  upon  the  theoretical  level  of  electronic  correlation. 


Table  5-8;  Relative  energies  of  minima  I and  II  and  of  transition  state  I-II  (Cs).  Eneigies  in 
kcal/mol. 


Theory  Level® 

Structure  I 

Structure  II 

Transition 

State 

Global 

Minimum 

UHF/3-21G 

8.68 

4.94 

9.11 

0.0 

UHF/6-31+G* 

9.78 

7.00 

11.58 

0.0” 

UMP2/6-31G* 

7.75 

9.00 

9.15 

4D 

O 

• 

o 

a:  AU  results  from  Ref.  [106]  except  UHF/3-21G,  b:  The  global  minimum  has  a Cj  structure. 

We  think  that  the  LiCNLi  molecule  is  a good  prototype  system  for  charge  transfer 
dynamics.  We  present  Mulliken  population  analysis,  dipole  moments,  and  isodensity 
surfaces  for  all  extrema  in  the  potential  energy  surface  of  LiCNLi.  Table  5-9  contains 
charges  and  dipole  moments  for  these  structures.  As  we  see  from  Table  5-9  almost  a 
complete  electron  transfer  takes  place  when  going  from  structure  I to  II  and  vice  versa. 
Since  the  atomic  charge  is  not  an  observable,  we  plot  in  Figures  5-41-5-44  the  electron 
isodensity  surfaces  of  these  structures  for  both  alpha  and  beta  spins,  and  confirm  what 
has  been  stated  using  ordinary  Mulliken  analysis.  We  should  mention  that  the  assignment 
of  alpha  or  beta  electrons  is  arbitrary,  and  we  choose  to  have  ten  alpha  electrons  and  9 
beta  electrons.  As  a result,  assuming  that  the  spin  contamination  is  small,  the  unpaired 
electron  is  alpha  and  this  explains  why  the  alpha  electron  density  is  asymmetric. 
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Table  5-9:  MuUiken  population  (a.u.)  and  dipx)le  moments  (D)  of  structures  I and  II,  transition 
state  I-II  and  global  minimum  at  UHF/3-21G  level. 


Structure  I 

Structure  II 

Transition 

State 

Global 

Minimum 

Li(C) 

0.704 

-0.168 

0.017 

0.283 

C 

0.061 

0.213 

0.138 

0.256 

N 

-0.599 

-0.761 

-0.674 

-0.821 

Li(N) 

-0.167 

0.716 

0.518 

0.283 

Dipole  (D) 

-16.74 

15.68 

-6.16 

-0.54 

alpha  density 


beta  density 


Figure  5-41:  Alpha  and  beta  isodensity  for  structure  I at  UHF/3-21G  level.  Isodensity 

is  0.003  a.u.  (1  a.u.  = 6.748  e/A3|.  Structure:  Li-C  = 197.98  pm  (=  1.980  A = 3.741 
a.u.),  C-N  = 115.02  pm  (=  1.150  A = 2.174  a.u.),  Li-N  = 192.59  pm  (=  1.926  A = 
3.639  a.u.).  SCF  energy:  -106.63313290  a.u.  = -2901.6368  eV 


alpha  density 


beta  density 


Figure  5-42:  Alpha  and  beta  isodensity  for  structure  II  at  UHF/3-21G  level.  Isodensity 

is  0.003  a.u.  (1  a.u.  = 6.748  e/A3|.  Structure:  Li-C  = 214.10  pm  (=  2.141  A = 4.046 
a.u.),  C-N  = 115.77  pm  (=  1.158  A = 2.188  a.u.),  Li-N  = 178.79  pm  (=  1.788  A = 
3.379  a.u.).  SCF  energy:  -106.63909258  a.u.  = -2901.7989  eV 


ti  cn 
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alpha  density 


beta  density 


Figure  5-43:  Alpha  and  beta  isodensity  of  the  transition  state  I-II  structure  (Cs)  at  UHF/3-21G 

level.  Isodensity  is  0.003  a.u.  (1  a.u.  = 6.748  e/A^).  Structure:  Li-C  = 208.75  pm  (=  2.088  A 
3.945  a.u.),  C-N  = 115.83  pm  (=  1.158  A = 2.189  a.u.),  Li-N  = 181.89  pm  (=  1.819  A = 
.437  a.u.),  Li-C-N  = 148.49°,  C-N-Li  = 138.61°.  SCF  energy:  -106.63244750  a.u.  = 
-2901.6180eV 


alpha  density 


beta  density 


Figure  5-44:  Alpha  and  beta  isodensity  of  the  global  minimum  structure  (C2v)  at  UHF/3-21G 

level.  Isodensity  is  0.003  a.u.  (1  a.u.  = 6.748  e/A^).  Structure:  C-N  = 117.46  pm  (=  1.175  A 
= 2.220  a.u.),  Li-N  = 189.22  pm  (=  1.892  A = 3.576  a.u.),  C-N-Li  = 137.12°.  SCF  energy: 
-106.64697039  a.u.  = -2902.0134  eV 
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The  choice  of  an  isodensity  value  of  0.003  a.u.  (1  a.u.  = I/oq  = 6.748  e/A^)  is 
based  on  the  fact  that  the  isodensity  in  the  range  of  0.002-0.003  a.u.  yields  molecular 
dimensions  that  agree  with  kinetic  theory  data,  so  the  expected  shape  of  the  molecule 
can  be  shown  [107].  Figures  5-41  and  5-42  corroborate  the  finding  of  the  Mulliken 
population  analysis  and  the  dipole  moment  calculations  for  structures  I and  II. 

We  are  unable  to  find  any  transition  state  that  is  linear.  Consequently,  the  reaction 
coordinate  for  the  electron  transfer  is  non-trivial  and  application  of  a two-state  type 
model  is  not  obvious.  As  a result,  we  decided  not  to  perform  any  dynamical  calculation 
and  instead  to  explore  the  intramolecular  electron  transfer  in  the  LiHLi  molecule.  The 
behavior  of  the  LiHLi  molecular  wave  function  and  potential  surface  is  common  among 
the  electron  transfer  systems  and  dynamical  results  from  the  END-SD-FGWP  formalism 
can  be  compared  with  other  theoretical  approaches. 

Structure,  Energetics,  and  Electronic  Population  of  LiHLi 

Before  presenting  the  results  for  LiHLi  we  should  specify  that  all  results  presented 
in  this  subsection  were  obtained  using  the  3-21G  basis  set  for  H [102]  and  3-214G  for 
Li  [108].  In  addition,  we  perform  time-independent  ab  initio  self-consistent  field  (SCF) 
calculations  to  determine  potential  surfaces  and  electronic  densities  of  the  LiHLi  and 
LiCNLi  molecules. 

The  linear  structure  of  the  open  shell  LiHLi  molecule  at  the  SCFAJHF  level  has 
a broken  symmetry  wave  function  which  is  strongly  (charge)  localized.  This  can  be 
seen  in  Figure  5-45,  where  we  plot  the  difference  of  Mulliken  charges  on  the  lithium 
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atoms  versus  the  relative  Li-H  bond  distances,  and  the  potential  energy  as  function  of 
the  reaction  coordinate  in  Figure  5-46. 


Figure  5-45:  Mulliken  charge  differences  at  SCFAJHF  level  in  the  Li(l)-H-Li(2)  molecule  as 
function  of  relative  bond  distances  (rj  - T2). 
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Figure  5-46:  Potential  energy  curves  of  the  Li(l)-H-Li(2)  molecule  at  SCF/UHF  level. 

As  we  can  see  from  Figures  5-45-5-46,  the  broken  symmetry  electronic  wave  function 
(lower  surface)  exhibits  an  abrupt  change  at  the  symmetric  structure  (ri  - r2  = 0,  where  ri 
= distance  Li(l)-H  and  r2  = distance  H-Li(2)),  This  leads  to  an  unphysical  discontinuous 
behavior  of  the  charge  density.  The  wave  function  describing  the  upper  surface  does  not 
show  this  behavior  perhaps  because  it  has  less  spin  contamination.  We  should  mention 
that  the  two  asymmetric  linear  structure  minima  are  in  fact  saddle  points,  since  they  have 
one  imaginary  vibrational  frequency  leading  to  the  bend  structure.  Nevertheless,  the  linear 
structures  are  a good  candidate  for  charge  transfer  studies  since  the  system  is  small  enough 
to  allow  us  to  perform  some  computational  experiments  within  the  ENDyne  program. 

Dynamics  of  Electron  Transfer  in  LiHLi 


We  present  the  dynamical  results  for  the  intramolecular  electron  transfer  in  LiHLi 
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obtained  using  the  ENDyne  program  within  the  END-SD-FGWP  formalism. 

The  procedure  utilized  here  is  the  following, 

(1)  make  a distortion  of  the  molecular  geometry  (contract  and/or  stretch  bonds)  relative 
to  the  linear  minima  structures; 

(2)  for  this  non-equilibrium  geometry  perform  a SCF/UHF  calculation  and  use  the 
converged  molecular  orbitals  as  the  initial  z parameters; 

(3)  let  the  system  evolve  in  time  and  analyze  the  molecular  properties  as  function  of  time. 

We  have  performed  time-dependent  calculations  for  four  different  initial  conditions 
(geometry  and  electronic  wave  function),  namely,  asymmetric  stretch  and  contraction  of 
the  Li-H  bonds  by  5.6%,  10%,  and  15%,  and  an  almost  symmetric  structure  (Iri  - T2\ 
= 0.004  a.u.).  Table  5-10  presents  the  structure  and  energies  of  these  initial  geometries 
for  evolution. 


Table  5-10:  MuUiken  population,  structure,  and  relative  energy  of  the  LiHLi  molecule  at 
SCF/UHF  level. 


Structure 

MuUiken  population 

Bond  distances  (a.u.) 

Relative 

energy 

(kcal/mol) 

Li(l) 

Li(2) 

Li(l)-H 

Li(2)-H 

Minimum 

2.387 

3.243 

3.0904 

3.4559 

0.0 

5.6% 

2.478 

3.195 

2.9170 

3.6491 

0.97 

10% 

2.498 

3.186 

2.7814 

3.8015 

3.2 

15% 

2.522 

3.178 

2.6110 

3.9743 

8.03 

SYM 

2.436 

3.219 

3.2549 

3.2586 

0.87 

The  results  of  the  time  evolution  of  these  structures  (Table  5-10)  are  presented  in 
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Figures  5-47-5-50.  The  changes  of  the  Li-H  bond  distances  as  functions  of  time  are 
plotted  in  Figure  5-47. 


Figure  5-47:  Bond  distances,  ri  = distance  Li(l)-H  and  tz  = distance  H-Li(2)  in  a.u.  as  a 
function  of  time.  ENDyne  calculation  with  a H/3-21G  and  Li/3-21+G  basis  set. 

Figure  5-48  has  the  results  for  the  Mulliken  population  of  the  lithium  atoms  during  the 
evolution,  and  the  differences  in  the  Mulliken  atomic  charges  are  plotted  in  Figure  5-49. 
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Figure  5-48:  Atomic  Mulliken  population  of  Li(l)  and  Li(2)  as  a function  of  time.  ENDyne 
calculation  with  a H/3-21G  and  Li/3-21+G  basis  set. 
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Figure  5-49:  Atomic  Mulliken  charge  differences  q[Li(l)]-q[Li(2)]  as  a function  of  time. 
ENDyne  calculation  with  a H/3-21G  and  Li/3-2 1+G  basis  set. 

Comparing  the  atomic  charges  as  functions  of  time  in  Figure  5-49  with  the  charges  in 
terms  of  the  reaction  coordinate  in  Figure  5-45  we  see  that  the  dynamical  treatment  seems 
to  better  describe  the  charge  transfer  since  it  does  not  exhibits  unphysical  discontinuity 
like  the  time-independent  calculation. 
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Figure  5-50:  Alpha  and  beta  spin  atomic  Mulliken  populations  of  Li(l)  and  Li(2)  in  the  15% 
stnicture  as  a function  of  time.  ENDyne  calculation  with  a H/3-21G  and  Li/3-2 1+G  basis  set. 

Figure  5-50  represents  a more  detailed  Mulliken  population  for  the  15%  structure. 
As  we  can  see  the  beta  spin  populations  remain  quite  constant  in  time  and  all  the  electron 
transfer  occur  between  the  alpha  spin  electronic  densities.  The  reason  for  this  is  that,  as 
we  said  before,  the  assignment  of  alpha  or  beta  spin  for  the  electrons  is  arbitrary,  and  as 
in  the  case  of  LiCNLi,  we  have  assigned  alpha  spin  to  the  unpaired  electron. 

Summary  and  Conclusion 


In  conclusion  of  this  chapter  we  again  point  out  that  the  application  of  the  END-SD- 
FGWP  method  to  electron  transfer  does  not  depend  upon  any  coupling  matrix  element 
neither  is  an  a priori  knowledge  of  the  reaction  coordinate  necessary.  In  addition  the 
ENDyne  results  for  intramolecular  charge  transfer  in  the  LiHLi  molecules  are  smooth  in 
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both  time  and  nuclear  displacement,  unlike  the  SCF/UHF  results  that  show  an  unphysical 
discontinuity  at  the  symmetric  structure. 


CHAPTER  6 
CONCLUSIONS 


The  exploration  of  the  capabilities  of  the  END-SD-FGWP  method  show  that  it  is 
successful  in  providing  an  appropriate  description  of  electron  nuclear  dynamics.  The 
main  advantages  of  the  END  formalism  are  that  it  (i)  does  not  need  the  knowledge  of 
potential  energy  surfaces,  (ii)  treats  appropriately  the  coupling  between  electrons  and 
nuclei,  and  (iii)  does  not  require  a priori  knowledge  of  the  reaction  coordinate  neither 
the  coupling  matrix  element  for  calculating  electron  transfer. 

The  approximations  employed,  namely  a single  determinantal  wave  function  for  the 
electrons  and  a classical  description  of  the  nuclei  seem  to  be  excellent  starting  points  for 
generating  qualitative  and  in  some  cases  even  quantitative  results  in  molecular  dynamics. 
For  instance,  in  the  case  of  the  intramolecular  electron  transfer  in  LiHLi,  the  Mulliken 
charges  generated  by  the  ab  initio  SCFAJHF  method  show  an  unphysical  discontinuous 
behavior  at  the  linear  symmetric  structure.  The  dynamical  calculations,  however,  do 
not  exhibit  this  behavior.  In  fact,  the  Mulliken  charges,  calculated  by  the  ENDyne 
program  based  upon  the  END-SG-FGWP  formalism,  change  with  a continuous  behavior. 
Although,  these  are  successful  results,  there  is  still  room  for  further  developments,  from 
both  formal  and  computational. 

Inclusion  of  electron  correlation  with  an  multireference  wave  function  [2]  would  bring 
the  END  formalism  to  a level  where,  predictive,  quantitative  results  could  be  obtained  for 
most  modest  size  molecular  system.  As  we  have  seen  from  the  scattering  results  of  the 
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+ He  and  + H2  systems  the  quantum  interferences  between  different  trajectories 
are  important  for  an  appropriate  description  of  the  differential  cross  sections.  We  present 
semiclassical  corrections  to  this  problem,  for  the  elastic  differential  cross  section,  which 
seem  to  work  well.  However,  for  inelastic  scattering  the  semiclassical  corrections  cannot 
be  straightforwardly  employed.  Thus,  an  appropriate  approach  within  the  END-SD- 
FGWP  formalism  needs  to  be  developed  in  order  to  take  into  account  these  quantum 
interference  effects.  The  proper  way  of  treating  this  problem  would  be  the  extension 
to  include  a quantum  description  of  the  nuclei.  In  other  words,  if  we  avoid  the  FGWP 
(classical)  approximation  we  could  solve  the  problem  of  quantum  interferences  in  the 
differential  cross  section  calculations.  This  quantum  treatment  of  the  nuclei  would  also 
provide  the  appropriate  approach  for  processes  such  as  nuclear  tunneling,  multi-channel 
reactions,  and  zero  point  energy  effects. 

The  ab  initio  implementation  of  the  END-SD-FGWP  method  requires  the  evaluation 
of  all  integrals  and  integral  derivatives  for  each  step  in  the  differential  equation  solver. 
This  makes  the  ab  initio  implementation  computationally  demanding,  and  limits  the 
applications  to  small  and  medium  sized  molecules.  We  have  presented  a simplification 
of  the  END-SD-FGWP  method  that  uses  the  NDDO  approximation  with  effective  cores. 
An  appropriate  parametrization  of  the  NDDO,  such  as  provided  by  the  AMI  and  PM3 
methods  has  proved  useful  and  accurate  in  describing  stationary  molecular  properties.  It 
is  our  expectation  that  the  same  quality  of  description  will  also  be  adequate  in  the  case  of 
time-dependent  methods.  From  a computational  point  of  view  the  savings  will  be  twofold: 
(i)  a decrease  by  several  orders  of  magnitude  of  the  number  of  molecular  integrals,  (ii) 


121 

the  elimination  of  high  frequency  motions  (core  electrons)  allowing  the  integrator  to  take 
larger  time  steps,  and  (iii)  a decrease  of  the  number  of  dynamical  equations. 
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